text:
C FINDBD.F - Gateway function for computing the initial state and/or
C            the matrices B and D of a discrete-time linear system,
C            given the system matrices A, C, and possibly B, D, and the
C            input and output trajectories, using SLICOT routine IB01CD.
C
C RELEASE 2.0 of SLICOT System Identification Toolbox.
C Based on SLICOT RELEASE 5.7. Copyright (c) 2002-2020 NICONET e.V.
C
C Matlab call:
C   [(x0)(,B(,D))(,V)(,rcnd)] = findBD(jobx0,comuse(,job),A(,B),C(,D),Y
C                                      (,U,tol,printw,ldwork))
C
C     [x0,B,V,rcnd] = findBD(1,1,1,A,C,Y,U)
C   [x0,B,D,V,rcnd] = findBD(1,1,2,A,C,Y,U)
C        [B,V,rcnd] = findBD(0,1,1,A,C,Y,U)
C      [B,D,V,rcnd] = findBD(0,1,2,A,C,Y,U)
C       [x0,V,rcnd] = findBD(1,2,1,A,B,C,Y,U)
C       [x0,V,rcnd] = findBD(1,2,2,A,B,C,D,Y,U)
C         [x0,rcnd] = findBD(0,2)                  (x0 = 0, rcnd = 1)
C       [x0,V,rcnd] = findBD(1,3,A,C,Y)
C
C Purpose:
C   To estimate the initial state and/or the matrices B and D of a
C   discrete-time system, given the system matrices A, C, and possibly
C   B, D, and the input and output trajectories of the system.
C
C Input parameters:
C   jobx0 - integer option to specify whether or not the initial state
C           should be computed:
C           = 1 : compute the initial state x(0);
C           = 0 : do not compute the initial state (possibly, because
C                 x(0) is known to be zero).
C   comuse- integer option to specify whether the system matrices B
C           and D should be computed or used:
C           = 1 : compute the matrices B and D, as specified by job;
C           = 2 : use the matrices B and D, as specified by job;
C           = 3 : do not compute/use the matrices B and D.
C   job   - integer option to determine which of the system matrices
C           B and D should be computed or used:
C           = 1 : compute/use the matrix B only (D is known to be zero);
C           = 2 : compute/use the matrices B and D.
C           job must not be specified if jobx0 = 0 and comuse = 2, or
C           if comuse = 3.
C   A     - the n-by-n system state matrix A.
C   B     - the n-by-m system input matrix B (input when jobx0 = 1 and
C           comuse = 2).
C   C     - the l-by-n system output matrix C.
C   D     - the l-by-m system matrix D (input when jobx0 = 1, comuse = 2
C           and job = 2).
C   Y     - the t-by-l output-data sequence matrix.  Column  j  of  Y
C           contains the  t  values of the j-th output component for
C           consecutive time increments.
C   U     - the t-by-m input-data sequence matrix (input when jobx0 = 1
C           and comuse = 2, or comuse = 1).  Column  j  of  U
C           contains the  t  values of the j-th input component for
C           consecutive time increments.
C   tol   - (optional) tolerance used for estimating the rank of
C           matrices. If  tol > 0,  then the given value of  tol  is
C           used as a lower bound for the reciprocal condition number;
C           an m-by-n matrix whose estimated condition number is less
C           than  1/tol  is considered to be of full rank.
C           Default: m*n*epsilon_machine where epsilon_machine is the
C           relative machine precision.
C   printw- (optional) switch for printing the warning messages.
C           = 1:  print warning messages;
C           = 0:  do not print warning messages.
C           Default:    printw = 0.
C   ldwork- (optional) the workspace size.
C           Default : computed by the formula
C           LDWORK = MAX( minimum workspace size needed, 2*CSIZE/3,
C                        CSIZE - ( m + l )*t - 2*n*( n + m + l ) - l*m )
C           where CSIZE is the cache size in double precision words.
C
C Output parameters:
C   x0    - the n-vector of estimated initial state x(0) (output when
C           jobx0 = 1, and set to 0 when jobx0 = 0 and comuse = 2).
C   B     - the n-by-m system input matrix B (output when comuse = 1).
C   D     - the l-by-m system matrix D (output when comuse = 1 and
C           job = 2).
C   V     - the n-by-n orthogonal matrix which reduces A to a real Schur
C           form (output when jobx0 = 1 or comuse = 1).
C   rcnd  - (optional) vector of length 1 (or 2, if comuse = 1 and
C           job = 2), containing the reciprocal condition numbers of
C           the matrices involved in rank decisions.
C
C Contributor:
C   V. Sima, Research Institute for Informatics, Bucharest, April 2000.
C
C Revisions:
C   V. Sima, July 2000, February 2001, September 2001, Apr. 2009,
C   Dec. 2012.
C
C **********************************************************************
C
      SUBROUTINE MEXFUNCTION( NLHS, PLHS, NRHS, PRHS )
C
C .. Parameters ..
      DOUBLE PRECISION  ZERO
      PARAMETER         ( ZERO = 0.0D0 )
C
C .. Mex-file interface parameters ..
      INTEGER           PLHS(*), PRHS(*)
      INTEGER*4         NLHS, NRHS
C
C .. Mex-file integer functions ..
      INTEGER           mxCreateDoubleMatrix, mxGetPr
      INTEGER*4         mxGetM, mxGetN, mxIsNumeric, mxIsComplex
C
C .. Scalar parameters used by SLICOT subroutines ..
      DOUBLE PRECISION  TOL
      INTEGER           INFO, IWARN, L, LDA, LDB, LDC, LDD, LDU, LDV,
     $                  LDWORK, LDY, M, N, NSMP
      CHARACTER         COMUSE, JOB, JOBX0
C
C .. Allocatable arrays ..
C !Fortran 90/95 (Fixed dimensions should be used with Fortran 77.)
      INTEGER,          ALLOCATABLE :: IWORK(:)
      DOUBLE PRECISION, ALLOCATABLE :: A(:,:), B(:,:), C(:,:), D(:,:),
     $                                 DWORK(:), U(:,:), V(:,:), X0(:),
     $                                 Y(:,:)
C
C .. Local variables and constant dimension arrays ..
      DOUBLE PRECISION  TEMP
      CHARACTER*120     TEXT
      INTEGER           CSIZE, CUSE, IC, IJOB, INI, IP, IPS, IQ, ISIZE,
     $                  ITMP, LDW1, LDW2, LDW3, LDWMIN, LIWORK, MA,
     $                  MINSMP, MINWLS, NA, NCOL, NCP1, TASK
      LOGICAL           PRINTW
C
C .. Cache size in double precision words ..
      PARAMETER          ( CSIZE = 64000 )
C
C .. External subroutines ..
      EXTERNAL          IB01CD
C
C ..Intrinsic functions..
      INTRINSIC         INT, MAX, MIN
C
C Check for proper number of arguments.
C
      IF ( NRHS.LT.2 ) THEN
         CALL mexErrMsgTxt
     $        ( 'findBD requires at least 2 input arguments' )
      ELSE IF ( NLHS.LT.1 ) THEN
         CALL mexErrMsgTxt
     $        ( 'findBD requires at least 1 output arguments' )
      END IF
C
C Check dimensions of input parameters and read/set scalar parameters.
C
C   jobx0
C
      IF ( mxGetM( PRHS(1) ).NE.1 .OR. mxGetN( PRHS(1) ).NE.1 ) THEN
         CALL mexErrMsgTxt( 'JOBX0 must be a scalar' )
      END IF
      IF ( mxIsNumeric( PRHS(1) ).EQ.0 .OR.
     $     mxIsComplex( PRHS(1) ).EQ.1 ) THEN
         CALL mexErrMsgTxt( 'JOBX0 must be an integer scalar' )
      END IF
      CALL mxCopyPtrToReal8( mxGetPr( PRHS(1) ), TEMP, 1 )
      TASK = TEMP
      IF ( TASK.LT.0 .OR. TASK.GT.1 ) THEN
         CALL mexErrMsgTxt
     $        ( 'JOBX0 has 0 or 1 the only admissible values' )
      END IF
C
      IF ( TASK.EQ.1 ) THEN
         JOBX0 = 'X'
      ELSE
         JOBX0 = 'N'
      END IF
C
C   comuse
C
      IF ( mxGetM( PRHS(2) ).NE.1 .OR. mxGetN( PRHS(2) ).NE.1 ) THEN
         CALL mexErrMsgTxt( 'COMUSE must be a scalar' )
      END IF
      IF ( mxIsNumeric( PRHS(2) ).EQ.0 .OR.
     $     mxIsComplex( PRHS(2) ).EQ.1 ) THEN
         CALL mexErrMsgTxt( 'COMUSE must be an integer scalar' )
      END IF
      CALL mxCopyPtrToReal8( mxGetPr( PRHS(2) ), TEMP, 1 )
      CUSE = TEMP
      IF ( CUSE.LT.1 .OR. CUSE.GT.3 ) THEN
         CALL mexErrMsgTxt
     $        ( 'COMUSE has 1, 2, or 3 the only admissible values' )
      END IF
C
      IF ( CUSE.EQ.1 ) THEN
         COMUSE = 'C'
      ELSE IF ( CUSE.EQ.2 ) THEN
         COMUSE = 'U'
      ELSE
         COMUSE = 'N'
      END IF
C
      IF ( TASK.EQ.1 .AND. CUSE.EQ.3 ) THEN
         IF ( NRHS.LT.5 ) THEN
            CALL mexErrMsgTxt
     $           ( 'findBD requires at least 5 input arguments' )
         END IF
      ELSE IF ( CUSE.EQ.1 ) THEN
         IF ( NRHS.LT.6 ) THEN
            CALL mexErrMsgTxt
     $           ( 'findBD requires at least 6 input arguments' )
         END IF
      END IF
C
C   job
C
      IF ( ( TASK.EQ.1 .AND. CUSE.EQ.2 ) .OR. CUSE.EQ.1 ) THEN
         IF ( mxGetM( PRHS(3) ).NE.1 .OR. mxGetN( PRHS(3) ).NE.1 ) THEN
            CALL mexErrMsgTxt( 'JOB must be a scalar' )
         END IF
         IF ( mxIsNumeric( PRHS(3) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(3) ).EQ.1 ) THEN
            CALL mexErrMsgTxt( 'JOB must be an integer scalar' )
         END IF
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(3) ), TEMP, 1 )
         IJOB = TEMP
         IF ( IJOB.LT.1 .OR. IJOB.GT.2 ) THEN
            CALL mexErrMsgTxt
     $           ( 'JOB has 1 or 2 the only admissible values' )
         END IF
         IP = 4
C
         IF ( IJOB.EQ.1 ) THEN
            JOB = 'B'
         ELSE
            JOB = 'D'
         END IF
      ELSE
         IP   = 3
         JOB  = 'B'
         IJOB = 1
      END IF
      IF ( TASK.EQ.1 .AND. CUSE.EQ.2 ) THEN
         IF ( IJOB.EQ.1 .AND. NRHS.LT.7 ) THEN
            CALL mexErrMsgTxt
     $           ( 'findBD requires at least 7 input arguments' )
         ELSE IF ( IJOB.EQ.2 .AND. NRHS.LT.8 ) THEN
            CALL mexErrMsgTxt
     $           ( 'findBD requires at least 8 input arguments' )
         END IF
      END IF
C
      IPS  = IP
      N    = 0
      M    = 0
      L    = 1
      NSMP = 1
C
C   A(n,n)
C
      IF ( NRHS.GE.IP ) THEN
         N  = mxGetM( PRHS(IP) )
         NA = mxGetN( PRHS(IP) )
         IF ( NA.NE.N ) THEN
            WRITE( TEXT, '(''A must have '',I6,'' rows and columns'')' )
     $             N
            CALL mexErrMsgTxt( TEXT )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
            CALL mexErrMsgTxt( 'A must be a real matrix' )
         END IF
         IP = IP + 1
      END IF
C
C   B(n,m)
C
      IF ( TASK.EQ.1 .AND. CUSE.EQ.2 .AND. NRHS.GE.IP ) THEN
         NA = mxGetM( PRHS(IP) )
         M  = mxGetN( PRHS(IP) )
         IF ( NA.NE.N ) THEN
            WRITE( TEXT, '(''B must have '',I6,'' rows'')' ) N
            CALL mexErrMsgTxt( TEXT )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
            CALL mexErrMsgTxt( 'B must be a real matrix' )
         END IF
         IP = IP + 1
      END IF
C
C   C(l,n)
C
      IF ( NRHS.GE.IP ) THEN
         L  = mxGetM( PRHS(IP) )
         NA = mxGetN( PRHS(IP) )
         IF ( NA.NE.N ) THEN
            WRITE( TEXT, '(''C must have '',I6,'' columns'')' ) N
            CALL mexErrMsgTxt( TEXT )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
            CALL mexErrMsgTxt( 'C must be a real matrix' )
         END IF
         IP = IP + 1
C
         IF ( L.LE.0 ) THEN
            CALL mexErrMsgTxt( 'The system has no outputs' )
         END IF
      END IF
C
C   D(l,m)
C
      IF ( TASK.EQ.1 .AND. CUSE.EQ.2 .AND. IJOB.EQ.2 .AND. NRHS.GE.IP )
     $      THEN
         NA = mxGetM( PRHS(IP) )
         MA = mxGetN( PRHS(IP) )
         IF ( NA.NE.L ) THEN
            WRITE( TEXT, '(''D must have '',I6,'' rows'')' ) L
            CALL mexErrMsgTxt( TEXT )
         END IF
         IF ( MA.NE.M ) THEN
            WRITE( TEXT, '(''D must have '',I6,'' columns'')' ) M
            CALL mexErrMsgTxt( TEXT )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
            CALL mexErrMsgTxt( 'D must be a real matrix' )
         END IF
         IP = IP + 1
      END IF
C
C  Y(txp), U(txm)
C
      IF ( NRHS.GE.IP ) THEN
         NSMP = mxGetM( PRHS(IP) )
         IF ( mxGetN( PRHS(IP) ).NE.L ) THEN
            WRITE( TEXT, '('' Y must have '',I6,'' columns'')' ) L
            CALL mexErrMsgTxt( TEXT )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
            CALL mexErrMsgTxt( 'Y must be a real matrix' )
         END IF
         IP = IP + 1
      END IF
C
      IF ( ( ( TASK.EQ.1 .AND. CUSE.EQ.2 ) .OR. CUSE.EQ.1 ) .AND.
     $       NRHS.GE.IP ) THEN
         MA = mxGetN( PRHS(IP) )
         IF ( CUSE.EQ.2 .AND. MA.NE.M ) THEN
            WRITE( TEXT, '(''U must have '',I6,'' columns'')' ) M
            CALL mexErrMsgTxt( TEXT )
         ELSE
            M = MA
         END IF
C
         IF ( M.GT.0 ) THEN
            IF ( mxGetM( PRHS(IP) ).NE.NSMP ) THEN
               CALL mexErrMsgTxt
     $             ( 'U must have the same row dimension as Y' )
            END IF
            IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $           mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
               CALL mexErrMsgTxt( 'U must be a real matrix' )
            END IF
         END IF
         IP = IP + 1
      END IF
C
      IF( CUSE.EQ.1 ) THEN
         NCOL = N*M
         IF( TASK.EQ.1 )
     $      NCOL = NCOL + N
         MINSMP = NCOL
         IF( IJOB.EQ.2 ) THEN
            MINSMP = MINSMP + M
            IQ     = MINSMP
         ELSE IF ( TASK.EQ.0 ) THEN
            IQ     = MINSMP
            MINSMP = MINSMP + 1
         ELSE
            IQ     = MINSMP
         END IF
      ELSE
         NCOL = N
         IF( TASK.EQ.1 ) THEN
            MINSMP = N
         ELSE
            MINSMP = 0
         END IF
         IQ = MINSMP
      END IF
C
      IF ( NSMP.LT.MINSMP ) THEN
         WRITE( TEXT, '('' The number of samples should be at least '',
     $                     I10)' ) MINSMP
         CALL mexErrMsgTxt( TEXT )
      END IF
C
C     tol
C
      TOL = ZERO
      IF ( NRHS.GE.IP ) THEN
         IF ( mxGetM( PRHS(IP) ).NE.1 .OR. mxGetN( PRHS(IP) ).NE.1 )
     $         THEN
            CALL mexErrMsgTxt( 'TOL must be a scalar' )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
              CALL mexErrMsgTxt( 'TOL must be a real scalar' )
         END IF
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), TOL, 1 )
         IP = IP + 1
      END IF
C
C     printw
C
      PRINTW = .FALSE.
      IF ( NRHS.GE.IP ) THEN
         IF ( mxGetM( PRHS(IP) ).NE.1 .OR.
     $        mxGetN( PRHS(IP) ).NE.1 ) THEN
            CALL mexErrMsgTxt( 'PRINTW must be a scalar' )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
            CALL mexErrMsgTxt( 'PRINTW must be an integer scalar' )
         END IF
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), TEMP, 1 )
         ITMP = TEMP
         IF ( ITMP.LT.0 .OR. ITMP.GT.1 ) THEN
            CALL mexErrMsgTxt
     $           ( 'PRINTW has 0 or 1 the only admissible values' )
         END IF
         PRINTW = ITMP.EQ.1
         IP = IP + 1
      END IF
C
C Determine the lenghts of working arrays.
C The default value for LDWORK is computed using the formula
C      LDWORK = MAX( minimum value needed, 2*CSIZE/3,
C     $              CSIZE - ( M + L )*NSMP - 2*N*( N + M + L ) - L*M )
C where CSIZE is the cache size in double precision words.
C If LDWORK is specified, but it is less than the minimum workspace size
C needed, that minimum value is used instead.
C
      LDA = MAX( 1, N )
      LDB = LDA
      LDC = L
      LDD = LDC
      LDV = LDA
      LDY = MAX( 1, NSMP )
      IF ( M.GT.0 .AND. ( ( TASK.EQ.1 .AND. CUSE.EQ.2 ) .OR.
     $                                      CUSE.EQ.1 ) ) THEN
         LDU = LDY
      ELSE
         LDU = 1
      END IF
C
      LIWORK = NCOL
      IF ( ( TASK.EQ.0 .AND. CUSE.NE.1 ) .OR. MAX( N, M ).EQ.0 ) THEN
         LDWMIN = 2
         LDWORK = 2
      ELSE
         IF ( IJOB.EQ.2 )
     $      LIWORK = MAX( LIWORK, M )
         IQ     = IQ*L
         NCP1   = NCOL + 1
         ISIZE  = NSMP*L*NCP1
         IF ( N.GT.0 .AND. TASK.EQ.1 ) THEN
            IF ( CUSE.EQ.1 ) THEN
               IC = 2*N*N + N
            ELSE
               IC = 2*N*N
            END IF
         ELSE
            IC = 0
         END IF
         MINWLS = NCOL*NCP1
         IF ( CUSE.EQ.1 ) THEN
            IF ( IJOB.EQ.2 )
     $         MINWLS = MINWLS + L*M*NCP1
            IF ( M.GT.0 .AND. IJOB.EQ.2 ) THEN
               IA = M + MAX( 2*NCOL, M )
            ELSE
               IA = 2*NCOL
            END IF
            LDW1 = N*N*M + MAX( IC, IA )
            IF ( TASK.EQ.1 )
     $         LDW1 = LDW1 + L*N
            LDW2 = ISIZE  + MAX( N + MAX( IC, IA ), 6*NCOL )
            LDW3 = MINWLS + MAX( IQ*NCP1 + LDW1, 6*NCOL )
            IF ( M.GT.0 .AND. IJOB.EQ.2 ) THEN
               LDW2 = MAX( LDW2, ISIZE  + 2*M*M + 6*M )
               LDW3 = MAX( LDW3, MINWLS + 2*M*M + 6*M )
               INI = 3
            ELSE
               INI = 2
            END IF
         ELSE
            ITAU = IC + L*N
            LDW2 = ISIZE  + 2*N + MAX( IC, 4*N )
            LDW3 = MINWLS + 2*N + MAX( IQ*NCP1 + ITAU, 4*N )
            INI  = 2
         END IF
         LDWMIN = INI + N*( N + M + L ) +
     $            MAX( 5*N, INI, MIN( LDW2, LDW3 ) )
         LDWORK = MAX( LDWMIN, 2*CSIZE/3, CSIZE - ( M + L )*NSMP -
     $                         2*N*( N + M + L ) - L*M )
      END IF
C
C   ldwork
C
      IF ( NRHS.GE.IP ) THEN
         IF ( mxGetM( PRHS(IP) ).NE.1 .OR.
     $        mxGetN( PRHS(IP) ).NE.1 ) THEN
            CALL mexErrMsgTxt( 'LDWORK must be a scalar' )
         END IF
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR.
     $        mxIsComplex( PRHS(IP) ).EQ.1 ) THEN
              CALL mexErrMsgTxt( 'LDWORK must be a real scalar' )
         END IF
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), TEMP, 1 )
         ITMP = INT( TEMP )
         IF ( ITMP.LT.LDWMIN ) THEN
            LDWORK = LDWMIN
         ELSE
            LDWORK = ITMP
         END IF
      END IF
C
C Allocate variable dimension local arrays.
C !Fortran 90/95
C
      ALLOCATE ( A( LDA, N ), B( LDB, M ), C( LDC, N ), D( LDD, M ),
     $           DWORK( LDWORK ), IWORK( LIWORK ), U( LDU, M ),
     $           V( LDV, N ), X0( N ), Y( LDY, L ) )
C
C Copy inputs from MATLAB workspace to locally allocated arrays.
C
      IP = IPS
      IF ( NRHS.GE.IP ) THEN
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), A, LDA*N )
         IP = IP + 1
         IF ( TASK.EQ.1 .AND. CUSE.EQ.2 ) THEN
            IF ( N.GT.0 ) THEN
               CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), B, LDB*M )
            END IF
            IP = IP + 1
         END IF
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), C, LDC*N )
         IP = IP + 1
         IF ( TASK.EQ.1 .AND. CUSE.EQ.2 .AND. IJOB.EQ.2 ) THEN
            CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), D, LDD*M )
            IP = IP + 1
         END IF
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), Y, NSMP*L )
         IP = IP + 1
         IF ( CUSE.EQ.1 .OR. ( TASK.EQ.1 .AND. CUSE.EQ.2 ) ) THEN
            CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), U, NSMP*M )
         END IF
      END IF
C
C Do the actual computations.
C
      CALL IB01CD( JOBX0, COMUSE, JOB, N, M, L, NSMP, A, LDA, B, LDB, C,
     $             LDC, D, LDD, U, LDU, Y, LDY, X0, V, LDV, TOL, IWORK,
     $             DWORK, LDWORK, IWARN, INFO )
      IF ( IWARN.NE.0 .AND. PRINTW ) THEN
         WRITE( TEXT, '(''  IWARN = '',I4,'' ON EXIT FROM IB01CD'')'
     $        ) IWARN
      END IF
      IF ( INFO.NE.0 ) THEN
         WRITE( TEXT, '(''INFO = '',I4,'' ON EXIT FROM IB01CD'')'
     $        ) INFO
      ELSE
C
C Copy output to MATLAB workspace.
C
         IP = 1
         IF ( TASK.EQ.1 .OR. ( TASK.EQ.0 .AND. CUSE.EQ.2 ) ) THEN
            PLHS(IP) = mxCreateDoubleMatrix( N, 1, 0 )
            CALL mxCopyReal8ToPtr( X0, mxGetPr( PLHS(IP) ), N )
            IP = IP + 1
         END IF
C
         IF ( CUSE.EQ.1 .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateDoubleMatrix( N, M, 0 )
            CALL mxCopyReal8ToPtr( B, mxGetPr( PLHS(IP) ), N*M )
            IP = IP + 1
            IF ( NLHS.GE.IP .AND. IJOB.EQ.2 ) THEN
               PLHS(IP) = mxCreateDoubleMatrix( L, M, 0 )
               CALL mxCopyReal8ToPtr( D, mxGetPr( PLHS(IP) ), L*M )
               IP = IP + 1
            END IF
         END IF
C
         IF ( ( TASK.EQ.1 .OR. CUSE.EQ.1 ) .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateDoubleMatrix( N, N, 0 )
            CALL mxCopyReal8ToPtr( V, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
C
         IF ( NLHS.GE.IP ) THEN
            NO = 1
            IF ( CUSE.EQ.1 .AND. M.GT.0 .AND. IJOB.EQ.2 )
     $         NO = 2
            PLHS(IP) = mxCreateDoubleMatrix( NO, 1, 0 )
            CALL mxCopyReal8ToPtr( DWORK(2), mxGetPr( PLHS(IP) ), NO )
         END IF
      END IF
C
C Deallocate local arrays.
C !Fortran 90/95
C
      DEALLOCATE( A, B, C, D, V, U, Y, IWORK, DWORK, X0 )
C
C Error and warning handling.
C
      IF ( IWARN.NE.0 .AND. PRINTW ) THEN
         CALL mexPrintf( TEXT )
      END IF
C
      IF ( INFO.NE.0 ) THEN
         CALL mexErrMsgTxt( TEXT )
      END IF
C
      RETURN
C *** Last line of findBD ***
      END

max_stars_repo_name: SLICOT/SLICOT-SystemIdentification

max_issues_repo_name: SLICOT/SLICOT-SystemIdentification

max_forks_repo_name: SLICOT/SLICOT-SystemIdentification


text:
      SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT)
C--
C--   SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT)
C--
C--   Purpose: compute saturation specific humidity and
C--            relative hum. from specific hum. (or viceversa)
C--   Input:   IMODE  : mode of operation
C--            NGP    : no. of grid-points
C--            TA     : abs. temperature
C--            PS     : normalized pressure   (=  p/1000_hPa) [if SIG < 0]
C--                   : normalized sfc. pres. (= ps/1000_hPa) [if SIG > 0]
C--            SIG    : sigma level
C--            QA     : specific humidity in g/kg [if IMODE > 0]
C--            RH     : relative humidity         [if IMODE < 0]
C--            QSAT   : saturation spec. hum. in g/kg
C--   Output:  RH     : relative humidity         [if IMODE > 0]
C--            QA     : specific humidity in g/kg [if IMODE < 0]
C--
      REAL TA(NGP), PS(*), QA(NGP), RH(NGP), QSAT(NGP)
C
C---  1. Compute Qsat (g/kg) from T (degK) and normalized pres. P (= p/1000_hPa)
C        If SIG > 0, P = Ps * sigma, otherwise P = Ps(1) = const.
C
      E0=  6.108E-3
      C1= 17.269
      C2= 21.875
      T0=273.16
      T1= 35.86
      T2=  7.66

      DO 110 J=1,NGP
        IF (TA(J).GE.T0) THEN
          QSAT(J)=E0*EXP(C1*(TA(J)-T0)/(TA(J)-T1))
        ELSE
          QSAT(J)=E0*EXP(C2*(TA(J)-T0)/(TA(J)-T2))
        ENDIF
  110 CONTINUE
C
      IF (SIG.LE.0.0) THEN
        DO 120 J=1,NGP
          QSAT(J)=622.*QSAT(J)/(PS(1)-0.378*QSAT(J))
  120   CONTINUE
      ELSE
        DO 130 J=1,NGP
          QSAT(J)=622.*QSAT(J)/(SIG*PS(J)-0.378*QSAT(J))
  130   CONTINUE
      ENDIF
C
C---  2. Compute rel.hum. RH=Q/Qsat (IMODE>0), or Q=RH*Qsat (IMODE<0)
C
      IF (IMODE.GT.0) THEN
        DO 210 J=1,NGP
          RH(J)=QA(J)/QSAT(J)
  210   CONTINUE
      ELSE IF (IMODE.LT.0) THEN
        DO 220 J=1,NGP
          QA(J)=RH(J)*QSAT(J)
  220   CONTINUE
      ENDIF
C
      RETURN
      END

      SUBROUTINE ZMEDDY (NLON,NLAT,FF,ZM,EDDY)
C
C *** Decompose a field into zonal-mean and eddy component
C
      REAL FF(NLON,NLAT), ZM(NLAT), EDDY(NLON,NLAT)
C
      RNLON=1./NLON
C
      DO 130 J=1,NLAT
C
        ZM(J)=0.
        DO 110 I=1,NLON
          ZM(J)=ZM(J)+FF(I,J)
 110    CONTINUE
        ZM(J)=ZM(J)*RNLON
C
        DO 120 I=1,NLON
          EDDY(I,J)=FF(I,J)-ZM(J)
 120    CONTINUE
C
 130  CONTINUE
C
C--
      RETURN
      END

ext: f


! { dg-do compile }
! { dg-options "" }
!
! PR fortran/32601
use, intrinsic :: iso_c_binding, only: c_loc, c_ptr
implicit none

! This was causing an ICE, but is an error because the argument to C_LOC
! needs to be a variable.
print *, c_loc(4) ! { dg-error "shall have either the POINTER or the TARGET attribute" }

end


text:
	PROGRAM TESTST
C************************************************************************
C* TESTST								*
C*									*
C* This program tests all of the subroutines in the string library.	*
C**									*
C* Log:									*
C* I. Graffman/RDS	 3/87	GEMPAK4					*
C* M. desJardins/GSFC	 5/89	Cleaned up ST_CRNM, ST_NUMB		*
C* S. Schotz/GSC	 1/90	Added ST_RMST				*
C* K. Brill/NMC		 9/91	Added ST_CKNR; ST_CKNI			*
C* M. desJardins/NMC	 8/92	Fixed input error to ST_ITOC		*
C* S. Jacobs/NMC	 8/94	Added ST_FLST				*
C* K. Tyle/GSC		 2/96	Increased character variable size	*
C* D. Plummer/NCEP	 3/96	Added ST_ILSR				*
C* S. Jacobs/NCEP	 3/96	Added ST_CLS2; ST_PELT			*
C* K. Tyle/GSC 		 7/96	Added ST_LSTC; ST_LSTI; ST_LSTF		*
C* D. Keiser/GSC	 8/96	Added ST_RXBL				*
C* S. Jacobs/NCEP	11/96	Added ST_RMNM				*
C* K. Tyle/GSC          12/96   Added ST_NULL				*
C* D.W.Plummer/NCEP      1/97   Added ST_RPST                           *
C* S. Maxwell/GSC	 4/97	Added ST_CRND				*
C* S. Jacobs/NCEP	 8/98	Added ST_CLSL				*
C* A. Hardy/GSC          8/98   Added ST_ETOA; ST_ATOE                  *
C* A. Hardy/GSC          9/98   Added a neg. number check for ST_ETOA   *
C* D.W.Plummer/NCEP	 3/00	Added ST_GTST				*
C* D. Kidwell/NCEP	 3/00	Added ST_RPSL				*
C* S. Jacobs/NCEP	 3/01	Added ST_SORT				*
C* D. Kidwell/NCEP	10/01	Added ST_WORD				*
C* M. Li/SAIC		02/02	Added ST_ISTP				*
C* T. Lee/SAIC		12/04	Added ST_OPCL				*
C************************************************************************
	CHARACTER	string*160, str1*160, strarr (10)*160, single*1
	CHARACTER	substr*160, str4 (10)*4, filarr (30)*72
        CHARACTER       repstr*160, outstr*160, instr*160, delim*1
	CHARACTER	first*10, last*10, inc*10, four (10)*4
	CHARACTER	inpsrt (20)*160, outsrt (20)*160
        BYTE            barray (60)
	INTEGER 	intg (10)
	REAL 		rlnum (10)
	LOGICAL 	abbr, istp
C------------------------------------------------------------------------
	CALL IN_BDTA ( iret )
	iostat = 0
	DO WHILE ( iostat .eq. 0 )
	string = ' '
	str1   = ' '
	WRITE (6, 10)
10	FORMAT (
     +  '  1 = ST_ABBR    2 = ST_ALNM    3 = ST_C2C     4 = ST_C2I'/
     +  '  5 = ST_C2R     6 = ST_CLST    7 = ST_CRNM    8 = ST_CTOI'/
     +  '  9 = ST_FIND   10 = ST_ILST   11 = ST_INCH   12 = ST_ITOC'/
     +  ' 13 = ST_ITOS   14 = ST_LCUC   15 = ST_LDSP   16 = ST_INLN'/
     +  ' 17 = ST_LSTR   18 = ST_NUMB   19 = ST_RANG   20 = ST_RLCH'/
     +  ' 21 = ST_RLST   22 = ST_RMBL   23 = ST_RNAN   24 = ST_STOI'/
     +  ' 25 = ST_UTAB   26 = ST_UNPR   27 = ST_FWRD   28 = ST_NXTS'/
     +  ' 29 = ST_INTG   30 = ST_RMST   31 = ST_FLST   32 = ST_ILSR'/
     +  ' 33 = ST_CLS2   34 = ST_PELT   35 = ST_LSTC   36 = ST_LSTI'/
     +  ' 37 = ST_LSTF   38 = ST_RXBL   39 = ST_RMNM   40 = ST_RPST'/
     +  ' 41 = ST_NULL   42 = ST_CRND   43 = ST_CLSL   44 = ST_ETOA'/
     +  ' 45 = ST_ATOE   46 = ST_GTST   47 = ST_RPSL   48 = ST_SORT'/
     +  ' 49 = ST_WORD   50 = ST_ISTP   51 = ST_OPCL')
	CALL TM_INT ( 'Select a subroutine number', .false.,
     +		       .false., 1, numsub, n, ier )
	IF ( ier .eq. 2 ) THEN
           iostat = -1
           numsub = -1
	END IF
2	FORMAT (A)
C------------------------------------------------------------------------
	IF (numsub .eq. 1) THEN
	    WRITE (6,*) 'Enter STRING: '
	    READ  (5,2)  string
	    WRITE (6,*) 'Enter STABBR:'
	    READ  (5,2)  str1
	    CALL ST_ABBR  ( string, str1, abbr, iret )
	    WRITE (6,*) 'ABBR, IRET: ', abbr, iret
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 2) THEN
	    WRITE (6,*) 'Enter CHRSTR:'
	    READ  (5,2)  single
	    CALL ST_ALNM  ( single, ityp, iret )
	    WRITE (6,*) 'ITYP, IRET: ', ityp, iret
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 3) THEN
	    WRITE (6,*) 'Enter STRING:'
	    READ  (5,2)  string
	    WRITE (6,*) 'Enter NEXP:'
	    READ  (5,*)  nexp
	    CALL ST_C2C  ( string, nexp, strarr, num, iret )
	    WRITE (6,*) 'IRET, NUM: ', iret, num
	    WRITE (6,*) ( strarr (i), ' ' ,i = 1, num)
	    CALL ER_WMSG  ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 4) THEN
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    WRITE (6,*)'Enter', nexp ,'integer strings with separators:'
	    READ (5,2) string
	    CALL ST_C2I (string, nexp, intg, num, iret)
	    WRITE (6,*) 'iret and number ret. ', iret, num
	    IF (num .GT.0) WRITE (6,*) (intg (i),i=1,num)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 5) THEN
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    WRITE (6,*)'Enter ', nexp ,'real strings with separators:'
	    READ (5, 2) string
	    CALL ST_C2R (string, nexp, rlnum, num, iret)
	    WRITE (6,*) 'iret, number returned ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (rlnum (i),i = 1,num)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 6) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6, *)' Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *)' Enter the default value:'
	    READ (5, 2) str1
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    CALL ST_CLST (string, single, str1, nexp, strarr,
     +	                  num, iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (strarr (i),i = 1,nexp)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 7) THEN
	    WRITE (6,*) 'Enter string:'
	    READ  (5,2)  string
	    CALL ST_CRNM  ( string, value, iret )
	    WRITE (6,*) 'VALUE, IRET: ', value, iret
	    CALL ER_WMSG  ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 8) THEN
	    WRITE (6, *)' Enter number of items (max 10):'
	    READ (5, *) nitem
	    WRITE (6, *)' Enter character array'
	    READ (5, 2) (four (i), i = 1, nitem)
	    CALL ST_CTOI (four, nitem, intg, iret)
	    WRITE (6, *) 'iret', iret
	    WRITE (6, *) (intg (i), i = 1, nitem)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 9) THEN
	    WRITE (6, *)' Enter number in list (max 10):'
	    READ (5, *) num
	    WRITE (6, *)' Enter string array:'
	    READ (5, 2) (strarr (i), i = 1, num)
	    WRITE (6, *) 'Enter string to find:'
	    READ (5, 2) string
	    CALL ST_FIND (string, strarr, num, ipos, iret)
	    WRITE (6, *)' Position = (0 = not found) ', ipos
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 10) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6, *)' Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *)' Enter the default value:'
	    READ (5, *) idef
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    CALL ST_ILST (string, single, idef, nexp, intg,
     +	                  num, iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (intg (i),i = 1,nexp)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 11) THEN
	    int=0
	    WRITE (6, *)'Enter an integer:'
	    READ (5, *) int
	    CALL ST_INCH (int, string, iret)
	    WRITE (6,*) 'string, iret ', string, iret
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 12) THEN
	    WRITE (6,*) 'Enter number of packed integers (max 10):'
	    READ  (5,*)  nexp
	    WRITE (6,*) 'Enter ', nexp, ' integers:'
	    READ  (5,*) (intg (i), i = 1, nexp)
	    CALL ST_ITOC  ( intg, nexp, str4, iret )
	    WRITE (6, 125) (str4 (i), i = 1, nexp)
125	    FORMAT (' ', 10(' ', A))
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 13) THEN
	    WRITE (6, *) 'Enter number of integers (max 10):'
	    READ (5,  *) nexp
	    WRITE (6, *)' Enter integer array:'
	    READ (5, *) (intg (i), i = 1, nexp)
	    CALL ST_ITOS (intg, nexp, nc, string, iret)
	    WRITE (6,*) ' iret and num chars: ', iret, nc
	    WRITE (6, *) ' string = ', string (1:nc)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 14) THEN
	    WRITE (6,*)' Enter a string: '
	    READ (5, 2) string
	    CALL ST_LCUC (string, str1, iret)
	    WRITE (6, *) str1
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 15) THEN
	    WRITE (6,*)'Enter a string'
	    READ (5,2) string
	    CALL ST_LDSP (string, str1, ncout, iret)
	    WRITE (6,*) 'ncout, string ', ncout
	    WRITE (6,*) str1
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 16) THEN
	    WRITE (6,*)'Enter an integer'
	    READ (5,*) int
	    CALL ST_INLN (int, string, lens, iret)
	    WRITE (6,*)' lens, string ', lens, string
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 17) THEN
	    WRITE (6,*) ' Enter a string: '
	    READ (5, 2) string
	    CALL ST_LSTR (string, isiz, iret)
            WRITE (6,*) 'isiz: ',isiz
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 18) THEN
	    WRITE (6,*) 'Enter STRING'
	    READ  (5,2)  string
	    CALL ST_NUMB ( string, number, iret )
	    WRITE (6 ,*) 'computed number, iret ', number, iret
	    CALL ER_WMSG  ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 19) THEN
	    WRITE (6, *) 'Enter a string range (use - ):'
	    READ (5, 2) string
	    CALL ST_RANG (string, first, last, inc, ityp, iret)
	    WRITE (6, *) ' first, last, inc, type, iret'
	    WRITE (6, *) first, last, inc, ityp, iret
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 20) THEN
	    WRITE (6,*)' Enter a real number and number of places:'
	    READ (5, *) rl, np
	    CALL ST_RLCH (rl, np, string, iret)
	    WRITE (6, *)' iret, string ', iret, string
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 21) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6, *)' Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *)' Enter the default value:'
	    READ (5, *) def
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    CALL ST_RLST (string, single, def, nexp, rlnum,
     +	                  num, iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (rlnum (i),i = 1,nexp)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 22) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    CALL ST_RMBL (string, str1, isiz, iret)
	    WRITE (6, *)' size and returned string: ', isiz
	    WRITE (6,*) str1
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 23) THEN
	    WRITE (6,*) 'Enter a string'
	    READ (5, 2) string
	    CALL ST_RNAN (string, str1, ncout, iret)
	    WRITE (6,*) 'string, nc, iret ', str1, ncout, iret
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 24) THEN
	    WRITE (6, *)' Enter string'
	    READ (5, 2) string
	    WRITE (6,*)' Enter number of characters in string:'
	    READ (5,*) nc
	    CALL ST_STOI (string, nc, nint, intg, iret)
	    WRITE (6, *) 'iret, nint ', iret, nint
	    IF (nint .gt. 0) WRITE (6, *) (intg (i), i = 1, nint)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 25) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5,2) string
	    WRITE (6,*) 'Enter number of characters in string'
	    READ (5,*)nc
	    CALL ST_UTAB (string, nc, str1, ier)
	    WRITE (6,*) str1
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 26) THEN
	    WRITE (6,*)' Enter string:'
	    READ (5, 2) string
	    WRITE (6,*) 'Enter length of string'
	    READ (5,*) nn
	    CALL ST_UNPR (string, nn, str1, lenout, iret)
	    WRITE (6,*)' lenout, outstr ', lenout, str1
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 27) THEN
	    WRITE (6, *) 'Enter a string:'
	    READ (5, 2) string
	    WRITE (6,*) 'Enter first, last, word number'
	    READ (5,*) ifirst, ilast, iword
	    CALL ST_FWRD (string, ifirst, ilast, iword, is, ie, iret)
	    WRITE (6,*)' start and end: ', is, ie
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 28) THEN
	    WRITE (6, *) 'Enter a string:'
	    READ (5, 2) string
	    WRITE (6,*) 'Enter first, last'
	    READ (5,*) ifirst, ilast
	    WRITE (6,*) 'Enter number of keywords (max = 10)'
	    READ (5, *) nexp
	    WRITE (6,*)' Enter ', nexp, 'keywords'
	    READ (5, 2) (strarr (i), i = 1, nexp)
	    DO i = 1, nexp
	        CALL ST_LSTR (strarr (i), intg (i), ier)
	    END DO
	    CALL ST_NXTS (string, ifirst, ilast, strarr, intg, nexp,
     +	                  ipos, istr, iret)
	    WRITE (6, *) ' ipos, istrg ', ipos, istr
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 29) THEN
	    WRITE (6, *)' Enter a string'
	    READ (5, 2) string
	    CALL ST_LSTR  ( string, lens, ier )
	    CALL ST_INTG (string (1:lens), num, iret)
	    WRITE (6,*) ' Integer = ', num
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 30) THEN
	    WRITE (6, *) ' Enter STRING'
	    READ  (5, 2)   string
            WRITE (6, *) ' Enter SUBSTR'
	    READ  (5, 2)   substr
	    CALL ST_RMST ( string, substr, ipos, string, iret )
	    WRITE (6,*)  'OUTSTR: ',string, ':'
	    WRITE (6,*)  ' IPOS, IRET: ', ipos, iret
	    IF (iret .ne. 0) CALL ER_WMSG ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 31) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6, *)' Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *)' Enter the default value:'
	    READ (5, 2) str1
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    CALL ST_FLST (string, single, str1, nexp, filarr,
     +	                  num, iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (filarr (i),i = 1,nexp)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 32) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6, *)' Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *)' Enter the default value:'
	    READ (5, *) idef
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    CALL ST_ILSR (string, single, idef, nexp, intg,
     +	                  num, iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (intg (i),i = 1,nexp)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 33) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6, *)' Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *)' Enter the default value:'
	    READ (5, 2) str1
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    CALL ST_CLS2 (string, single, str1, nexp, strarr,
     +	                  num, iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (strarr (i),i = 1,nexp)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 34) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6,*)' Enter item number:'
	    READ (5,*) nelmt
	    WRITE (6, *)' Enter the default value:'
	    READ (5, 2) str1
	    CALL ST_PELT (string, nelmt, str1, num, substr,
     +	                  iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    WRITE (6,*) substr
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
	  ELSE IF ( numsub .eq. 35 ) THEN
	    WRITE (6, *) 'Enter number of elements in the array:'
	    READ (5, *) num
	    DO I = 1, num
	      WRITE (6, *) 'Enter string #',i,':'
	      READ (5, 2) strarr ( i )
	    END DO
	    WRITE (6, *) 'Enter a separator:'
	    READ (5, 2) single
	    CALL ST_LSTC ( strarr, num, single, string, iret )
	    WRITE ( 6, * ) 'String: ', string
	    WRITE ( 6, * ) 'iret = ', iret
	    IF ( iret .ne. 0 ) CALL ER_WMSG  ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF ( numsub .eq. 36 ) THEN
	    WRITE (6, *) 'Enter number of elements in the array:'
	    READ (5, *) num
	    DO I = 1, num
	      WRITE (6, *) 'Enter integer #',i,':'
	      READ (5, *) intg ( i )
	    END DO
	    WRITE (6, *) 'Enter a separator:'
	    READ (5, 2) single
	    CALL ST_LSTI ( intg, num, single, string, iret )
	    WRITE ( 6, * ) 'String: ', string
	    WRITE ( 6, * ) 'iret = ', iret
	    IF ( iret .ne. 0 ) CALL ER_WMSG  ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF ( numsub .eq. 37 ) THEN
	    WRITE (6, *) 'Enter number of elements in the array:'
	    READ (5, *) num
	    DO I = 1, num
	      WRITE (6, *) 'Enter real #',i,':'
	      READ (5, *) rlnum ( i )
	    END DO
	    WRITE (6, *) 'Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *) 'Enter number of decimal places to include:'
	    READ (5, *) np
	    CALL ST_LSTF ( rlnum, num, single, np, string, iret )
	    WRITE ( 6, * ) 'String: ', string
	    WRITE ( 6, * ) 'iret = ', iret
	    IF ( iret .ne. 0 ) CALL ER_WMSG  ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 38) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    CALL ST_RXBL (string, str1, isiz, iret)
	    WRITE (6, *)' size and returned string: ', isiz
	    WRITE (6,*) str1
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 39) THEN
	    WRITE (6,*) 'Enter a string'
	    READ (5, 2) string
	    CALL ST_RMNM (string, str1, nncr, ncout, iret)
	    WRITE (6,*) 'string, nncr, nc, iret: ',
     +			str1, nncr, ncout, iret
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 40) THEN
            WRITE (6, *) ' Enter STRING'
            READ  (5, 2)   string
            WRITE (6, *) ' Enter SUBSTR'
            READ  (5, 2)   substr
            WRITE (6, *) ' Enter REPSTR'
            READ  (5, 2)   repstr
            CALL ST_RPST ( string, substr, repstr, ipos, string, iret )
            WRITE (6,*)  'OUTSTR: ',string
            WRITE (6,*)  ' IPOS, IRET: ', ipos, iret
            IF (iret .ne. 0) CALL ER_WMSG ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 41) THEN
            WRITE (6,*) 'Enter a string'
            READ (5, 2) string
            CALL ST_NULL (string, str1, isiz, iret)
            WRITE (6, *)' size and returned string: ', isiz
            WRITE (6,*) str1
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 42) THEN
            WRITE (6,*) 'Enter string:'
            READ  (5,2)  string
            CALL ST_CRND  ( string, value, ndec, iret )
            WRITE (6,*) 'VALUE, NDEC, IRET: ', value, ndec, iret
            CALL ER_WMSG  ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 43) THEN
	    WRITE (6, *)' Enter a string:'
	    READ (5, 2) string
	    WRITE (6, *)' Enter a separator:'
	    READ (5, 2) single
	    WRITE (6, *)' Enter the default value:'
	    READ (5, 2) str1
	    WRITE (6,*)' Enter number expected'
	    READ (5,*) nexp
	    CALL ST_CLSL (string, single, str1, nexp, strarr,
     +	                  num, iret)
	    WRITE (6,*) 'iret, number in list ', iret, num
	    IF (num .gt. 0) WRITE (6,*) (strarr (i),i = 1,nexp)
	    IF (iret .ne. 0) CALL ER_WMSG ('ST', iret, ' ', ier)
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 44) THEN
	    WRITE (6, *) 'Enter the number of elements in the array:'
	    READ (5, *) num
	    DO I = 1, num
	        WRITE (6, *)' Enter EBCDIC code # ',i,':'
	        READ (5, *) barray(I)
	    END DO
            CALL ST_ETOA ( barray, num, outstr, iret )
            IF (num .gt. 0)  WRITE(6,*)'ASCII string:  ',outstr(1:num)
            IF ( num .lt. 0 ) WRITE(6,*)'Invalid number of elements.'
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 45) THEN
	    WRITE (6, *)' Enter the ASCII characters.'
	    READ (5, 2) instr
            CALL ST_LSTR ( instr, num, iret )
            CALL ST_ATOE ( instr, num, barray, iret )
            write(6,*)'NUM ', num
            write(6,*)'BYTE values: ', (barray(i), i = 1, num )
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 46) THEN
            WRITE (6, *) ' Enter STRING (template)'
            READ  (5, 2)   string
            WRITE (6, *) ' Enter SUBSTR'
            READ  (5, 2)   substr
            WRITE (6, *) ' Enter DELIMITER'
            READ  (5, 2)   delim
            WRITE (6, *) ' Enter STRING (extract)'
            READ  (5, 2)   repstr
            CALL ST_GTST ( string, substr, delim, repstr,
     +			   string, len, iret )
            WRITE (6,*)  'OUTSTR: ',string
            WRITE (6,*)  'LENGTH: ',len
            WRITE (6,*)  ' IRET: ', iret
            IF (iret .ne. 0) CALL ER_WMSG ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 47) THEN
            WRITE (6, *) ' Enter STRING'
            READ  (5, 2)   string
            WRITE (6, *) ' Enter SUBSTR'
            READ  (5, 2)   substr
            WRITE (6, *) ' Enter LENSUB'
            READ  (5, *)   lensub
            WRITE (6, *) ' Enter REPSTR'
            READ  (5, 2)   repstr
            WRITE (6, *) ' Enter LENREP'
            READ  (5, *)   lenrep
            CALL ST_RPSL ( string, substr, lensub, repstr, lenrep,
     +			   ipos, string, iret )
            WRITE (6,*)  'OUTSTR: ',string
            WRITE (6,*)  ' IPOS, IRET: ', ipos, iret
            IF (iret .ne. 0) CALL ER_WMSG ( 'ST', iret, ' ', ier )
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 48) THEN
            WRITE (6, *) ' Enter the type of sort'
            READ  (5, *)   itype
            WRITE (6, *) ' Enter number of input strings (<=20)'
            READ  (5, *)   nstr
	    DO  i = 1, nstr
		WRITE (6, *) ' Enter STRING', i
		READ  (5, 2)   inpsrt(i)
	    END DO
	    CALL ST_SORT ( itype, nstr, inpsrt, nout, outsrt, iret )
            WRITE (6,*)  ' IRET: ', iret
            IF (iret .ne. 0) CALL ER_WMSG ( 'ST', iret, ' ', ier )
	    WRITE (6,*)  ' Sorted array:'
	    DO  i = 1, nout
	    	WRITE (6,*) outsrt(i)
	    END DO
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 49) THEN
	    WRITE (6,*) 'Enter CHRSTR:'
	    READ  (5,2)  string
	    CALL ST_WORD  ( string, ityp, iret )
	    WRITE (6,*) 'ITYP, IRET: ', ityp, iret
C------------------------------------------------------------------------
          ELSE IF (numsub .eq. 50) THEN
            WRITE (6,*) 'Enter STRING: '
            READ  (5,2)  string
            CALL ST_ISTP  ( string, istp, iret )
            WRITE (6,*) 'ISTP, IRET: ', istp, iret
C------------------------------------------------------------------------
	  ELSE IF (numsub .eq. 51) THEN
	    WRITE (6, *) ' Enter INSTR '
            READ  (5, 2)   instr
            WRITE (6, *) ' Enter LOCOPN '
            READ  (5, *)   locopn
            CALL ST_OPCL ( instr, locopn, outstr, loccls, iret )
            WRITE (6,*)  'OUTSTR: ', outstr
            WRITE (6,*)  'LOCCLS: ', loccls
            WRITE (6,*)  ' IRET: ', iret
            IF (iret .ne. 0) CALL ER_WMSG ( 'ST', iret, ' ', ier )
	END IF
	END DO
C*
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C$Procedure PCKOPN ( PCK, open new file )

      SUBROUTINE PCKOPN ( NAME, IFNAME, NCOMCH, HANDLE )

C$ Abstract
C
C     Create a new PCK file, returning the handle of the opened file.
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     PCK
C
C$ Keywords
C
C     PCK
C
C$ Declarations

      IMPLICIT NONE

      CHARACTER*(*)         NAME
      CHARACTER*(*)         IFNAME
      INTEGER               NCOMCH
      INTEGER               HANDLE

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     NAME       I   The name of the PCK file to be opened.
C     IFNAME     I   The internal filename for the PCK.
C     NCOMCH     I   The number of characters to reserve for comments.
C     HANDLE     O   The handle of the opened PCK file.
C
C$ Detailed_Input
C
C     NAME     The name of the PCK file to be created.
C
C     IFNAME   The internal filename for the PCK file that is being
C              created. The internal filename may be up to 60 characters
C              long. If you do not have any conventions for tagging your
C              files, an internal filename of 'PCK_file' is perfectly
C              acceptable. You may also leave it blank if you like.
C
C     NCOMCH   This is the space, measured in characters, to be
C              initially set aside for the comment area when a new PCK
C              file is opened. The amount of space actually set aside
C              may be greater than the amount requested, due to the
C              manner in which comment records are allocated in an PCK
C              file. However, the amount of space set aside for comments
C              will always be at least the amount that was requested.
C
C              The value of NCOMCH should be greater than or equal to
C              zero, i.e., 0 <= NCOMCH. A negative value, should one
C              occur, will be assumed to be zero.
C
C$ Detailed_Output
C
C     HANDLE   The handle of the opened PCK file. If an error occurs
C              when opening the file, the value of this variable should
C              not be used, as it will not represent a valid handle.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If the value of NCOMCH is negative, a value of zero (0) will
C        be used for the number of comment characters to be set aside
C        for comments.
C
C     2) If an error occurs while attempting to open a CK file the
C        value of HANDLE will not represent a valid file handle.
C
C$ Files
C
C     See NAME and HANDLE.
C
C$ Particulars
C
C     Open a new PCK file, reserving room for comments if requested.
C
C$ Examples
C
C     Suppose that you want to create a new PCK file called 'new.PCK'
C     that contains a single type 2 PCK segment and has room for at
C     least 5000 comment characters. The following code fragment should
C     take care of this for you, assuming that all of the variables
C     passed to the PCK type 2 segment writer have appropriate values.
C
C        NAME   = 'new.pck'
C        IFNAME = 'Test PCK file'
C
C        CALL PCKOPN ( NAME, IFNAME, 5000, HANDLE )
C        CALL PCKW02 ( HANDLE, BODY,   FRAME,  FIRST,   LAST,
C     .                SEGID,  INTLEN, N,      POLYDG,  CDATA,
C     .                BTIME                                   )
C        CALL PCKCLS ( HANDLE )
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C      None.
C
C$ Author_and_Institution
C
C      K.R. Gehringer    (JPL)
C
C$ Version
C
C-    SPICELIB Version 2.0.1, 03-JAN-2014 (EDW)
C
C        Minor edits to Procedure; clean trailing whitespace.
C        Corrected order of header sections to conform to NAIF
C        standard.
C
C-    SPICELIB Version 2.0.0, 09-NOV-2006 (NJB)
C
C        Routine has been upgraded to support comment
C        area allocation using NCOMCH.
C
C-    SPICELIB Version 1.0.0, 26-JAN-1995 (KRG)
C
C-&

C$ Index_Entries
C
C     open a new pck file
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               FAILED
      LOGICAL               RETURN
C
C     Local parameters
C
      CHARACTER*(*)         TYPE
      PARAMETER           ( TYPE = 'PCK'    )

C
C     DAF ND and NI values for PCK files.
C
      INTEGER               ND
      PARAMETER           ( ND = 2 )

      INTEGER               NI
      PARAMETER           ( NI = 5 )
C
C     Length of a DAF comment record, in characters.
C
      INTEGER               MXCREC
      PARAMETER           ( MXCREC = 1000 )

C
C     Local variables
C
      INTEGER               NCOMR

C
C     Standard SPICE error handling.
C
      IF ( RETURN() ) THEN
         RETURN
      END IF

      CALL CHKIN ( 'PCKOPN' )
C
C     Compute the number of comment records that we want to allocate, if
C     the number of comment characters requested is greater than zero,
C     we always allocate an extra record to account for the end of line
C     marks in the comment area.
C
      IF ( NCOMCH .GT. 0 ) THEN
         NCOMR = ( NCOMCH - 1 )/MXCREC  +  1
      ELSE
         NCOMR = 0
      END IF

C
C     Just do it. All of the error handling is taken care of for us.
C
      CALL DAFONW( NAME, TYPE, ND, NI, IFNAME, NCOMR, HANDLE )

      IF ( FAILED() ) THEN
C
C        If we failed, make sure that HANDLE does not contain a value
C        that represents a valid DAF file handle.
C
         HANDLE = 0

      END IF

      CALL CHKOUT ( 'PCKOPN' )
      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    @(#)ldptizon.f	20.9 3/29/99
C****************************************************************
C
C     File: ldptizon.f
C
C     Purpose: Routine to load PTI LTC Zone data from raw data file
C
c     Return code:  n = 0 : Success
c                   n = 1 : Error
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_pti.f
C
C****************************************************************
      integer function ldptizon (xbuf, options, numver, error)
      integer error, options(*), numver
      character *(*) xbuf

      include 'ipfinc/parametr.inc'

      include 'ipfinc/prt.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/arcntl.inc'

      character word(20)*10, zonename*2
      logical first, finished
      integer ftn_atoi, ptihashz, fnd_ptiy, status

      ldptizon = 0
      error = 0
      last = index (xbuf, '/') - 1
      if (last .le. 0) last = lastch(xbuf)
      call ptiscan (xbuf(1:last), word, nwrd, ',',  ' ')

      do i = nwrd+1, 2
        word(i) = ' '
      enddo
      npti1 = ftn_atoi(word(1))
      if (npti1 .eq. 0) then
        ldptizon = 1
      else
c
c       Generate unique zonename if necessary
c
        zonename = word(2)
        status = fnd_ptiy (zonename)
        finished = (status .eq. 0)
        next = 48
        do while (next .le. 90 .and. .not. finished)
           zonename = zonename(1:1) // char(next)
           status = fnd_ptiy (zonename)
           if (status .eq. 0) then
              finished = .true.
           else if (next .eq. 57) then
              next = 35
           else if (next .ge. 35 .and. next .le. 37) then
              next = next + 1
           else if (next .eq. 38) then
              next = 42
           else if (next .eq. 42) then
              next = 43
           else if (next .eq. 43) then
              next = 45
           else if (next .eq. 45) then
              next = 46
           else if (next .eq. 46) then
              next = 65
           else
              next = next + 1
           endif
        enddo
        if (zonename .ne. word(2)(1:2)) then
          write (errbuf(1), 10000) npti1, word(2), zonename
10000     format (' Duplicate zone name (1:2) on PTI Zone record (',
     &       i6, 1x, a, ') renamed to (', a, ')')
          call prterx ('W',1)
        endif

        num = ptihashz (npti1, zonename)
        if (num .le. 0) then
          write (errbuf(1), 10010) npti1, zonename
10010     format (' Zone on PTI Zone record (',
     &       i6, 1x, a, ') could not be hashed.')
          call prterx ('W',1)
          error = 1
          go to 900
        endif
      endif
  900 continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


module advance_module

  use bl_error_module
  use kernels_module
  use multifab_module
  use time_module
  use transport_properties
  use variables_module

  implicit none

  double precision, save, public :: wt_fillboundary = 0.d0
  double precision, save, public :: wt_ctoprim = 0.d0
  double precision, save, public :: wt_chemterm = 0.d0
  double precision, save, public :: wt_transprop = 0.d0
  double precision, save, public :: wt_hypdiff = 0.d0

  private
  public advance

contains

  subroutine advance(U, dt, courno, dx, istep)
    use smcdata_module
    type(multifab),    intent(inout) :: U
    double precision,  intent(inout) :: dt, courno
    double precision,  intent(in   ) :: dx(3)
    integer, intent(in) :: istep

    call multifab_setval(Unew, 0.d0, .true.)

    call dUdt(U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
    call set_dt(dt, courno, istep)
    call update_rk3(Zero,Unew, One,U, dt,Uprime)
    call reset_density(Unew)

    call dUdt(Unew, Uprime, Q, mu, xi, lam, Ddiag, dx)
    call update_rk3(OneQuarter, Unew, ThreeQuarters, U, OneQuarter*dt, Uprime)
    call reset_density(Unew)

    call dUdt(Unew, Uprime, Q, mu, xi, lam, Ddiag, dx)
    call update_rk3(OneThird, U, TwoThirds, Unew, TwoThirds*dt, Uprime)
    call reset_density(U)
  end subroutine advance

  !
  ! Compute new time-step size
  !
  subroutine set_dt(dt, courno, istep)

    use probin_module, only : cflfac, fixed_dt, init_shrink, max_dt, max_dt_growth, small_dt, stop_time

    double precision, intent(inout) :: dt
    double precision, intent(in   ) :: courno
    integer,          intent(in   ) :: istep

    double precision :: dtold

    if (fixed_dt > 0.d0) then

       dt = fixed_dt

       if (parallel_IOProcessor()) then
          print*, ""
          print*, "   Setting fixed dt =",dt
          print*, ""
       end if

    else

       dtold = dt
       dt    = cflfac / courno

       if (parallel_IOProcessor()) then
          print*, "   CFL: dt =", dt
       end if

       if (istep .eq. 1) then
          dt = dt * init_shrink
          if (parallel_IOProcessor()) then
             print*,'   Limited by init_shrink: dt =',dt
          end if
       else
          if (dt .gt. dtold * max_dt_growth) then
             dt = dtold * max_dt_growth
             if (parallel_IOProcessor()) then
                print*,'   Limited by dt_growth: dt =',dt
             end if
          end if
       end if

       if(dt .gt. max_dt) then
          if (parallel_IOProcessor()) then
             print*,'   Limited by max_dt: dt =',max_dt
          end if
          dt = max_dt
       end if

       if (dt < small_dt) then
          call bl_error("ERROR: timestep < small_dt")
       endif

       if (stop_time > 0.d0) then
          if (time + dt > stop_time) then
             dt = stop_time - time
             if (parallel_IOProcessor()) then
                print*, "   Limited by stop_time: dt =",dt
             end if
          end if
       end if

    end if

  end subroutine set_dt

  !
  ! Compute U1 = a U1 + b U2 + c Uprime.
  !
  subroutine update_rk3 (a,U1,b,U2,c,Uprime)

    type(multifab),   intent(in   ) :: U2, Uprime
    type(multifab),   intent(inout) :: U1
    double precision, intent(in   ) :: a, b, c

    integer :: lo(3), hi(3), i, j, k, m, n, nc
    type(mfiter) :: mfi
    type(box) :: tbx
    double precision, pointer, dimension(:,:,:,:) :: u1p, u2p, upp

    nc = ncomp(U1)

    !$omp parallel private(i,j,k,m,n,lo,hi,mfi,tbx,u1p,u2p,upp)

    call mfiter_build(mfi, U1, tiling=.true.)

    do while (next_tile(mfi,n))

       tbx = get_tilebox(mfi)
       lo = lwb(tbx)
       hi = upb(tbx)

       u1p => dataptr(U1,    n)
       u2p => dataptr(U2,    n)
       upp => dataptr(Uprime,n)

       do m = 1, nc
          do k = lo(3),hi(3)
             do j = lo(2),hi(2)
                do i = lo(1),hi(1)
                   u1p(i,j,k,m) = a * u1p(i,j,k,m) + b * u2p(i,j,k,m) + c * upp(i,j,k,m)
                end do
             end do
          end do
       end do
    end do
    !$omp end parallel

  end subroutine update_rk3

  !
  ! Compute dU/dt given U.
  !
  ! The Courant number (courno) is also computed if passed.
  !
  subroutine dUdt (U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)

    use probin_module, only: cfl_int, fixed_dt, do_tiling

    type(multifab),   intent(inout) :: U, Uprime, Q, mu, xi, lam, Ddiag
    double precision, intent(in   ) :: dx(3)
    integer,          intent(in   ), optional :: istep
    double precision, intent(inout), optional :: courno

    integer :: lo(3), hi(3)
    integer :: n, ng
    type(mfiter) :: mfi
    type(box) :: tbx

    logical :: update_courno
    double precision :: courno_proc

    integer :: qlo(4), qhi(4), uplo(4), uphi(4), ulo(4), uhi(4)
    double precision, pointer, dimension(:,:,:,:) :: up, qp, mup, xip, lamp, Ddp, upp

    double precision :: wt1, wt2

    update_courno = .false.
    if (present(courno) .and. present(istep) .and. fixed_dt.le.0.d0) then
       if (mod(istep,cfl_int).eq.1 .or. cfl_int.le.1) then
          update_courno = .true.
       end if
    end if

    wt1 = parallel_wtime()
    call multifab_fill_boundary(U)
    wt_fillboundary = wt_fillboundary + (parallel_wtime()-wt1)

    call multifab_setval(Uprime, 0.d0)

    !
    ! Calculate primitive variables based on U
    !
    wt1 = parallel_wtime()
    call ctoprim(U, Q)
    wt_ctoprim = wt_ctoprim + (parallel_wtime()-wt1)

    if (update_courno) then
       courno_proc = -1.d50
       call compute_courno(Q, dx, courno_proc)
    end if

    !
    ! chemistry
    !
    wt1 = parallel_wtime()
    do n=1,nfabs(Q)

       qp  => dataptr(Q,n)
       upp => dataptr(Uprime,n)

       qlo = lbound(qp)
       qhi = ubound(qp)
       uplo = lbound(upp)
       uphi = ubound(upp)

       lo = lwb(get_box(Q,n))
       hi = upb(get_box(Q,n))

       call chemterm_3d(lo,hi,qp,qlo(1:3),qhi(1:3),upp,uplo(1:3),uphi(1:3))
    end do
    wt2 = parallel_wtime()
    wt_chemterm = wt_chemterm + (wt2-wt1)

    !
    ! transport coefficients
    !
    call get_transport_properties(Q, mu, xi, lam, Ddiag)
    wt1 = parallel_wtime()
    wt_transprop = wt_transprop + (wt1-wt2)

    !
    ! Hyperbolic and Transport terms
    !
    if (do_tiling) then
       !$omp parallel private(n,lo,hi,up,ulo,uhi,upp,uplo,uphi) &
       !$omp private(qp,qlo,qhi,mup,xip,lamp,Ddp,mfi,tbx)

       call mfiter_build(mfi, Q, tiling=.true.)

       do while (next_tile(mfi,n))

          tbx = get_tilebox(mfi)
          lo = lwb(tbx)
          hi = upb(tbx)

          up => dataptr(U,n)
          upp=> dataptr(Uprime,n)
          qp => dataptr(Q,n)
          mup  => dataptr(mu   , n)
          xip  => dataptr(xi   , n)
          lamp => dataptr(lam  , n)
          Ddp  => dataptr(Ddiag, n)

          ulo = lbound(up)
          uhi = ubound(up)
          qlo = lbound(qp)
          qhi = ubound(qp)
          uplo = lbound(upp)
          uphi = ubound(upp)

          call narrow_diffterm_3d(lo,hi,dx,qp,qlo(1:3),qhi(1:3),upp,uplo(1:3),uphi(1:3), &
               mup,xip,lamp,Ddp)

          call hypterm_3d(lo,hi,dx,up,ulo(1:3),uhi(1:3),qp,qlo(1:3),qhi(1:3),&
               upp,uplo(1:3),uphi(1:3))

       end do
       !$omp end parallel
    else

       ng = nghost(U)

       do n=1,nfabs(Q)
          qp  => dataptr(Q,n)
          up  => dataptr(U,n)
          upp => dataptr(Uprime,n)

          mup  => dataptr(mu   , n)
          xip  => dataptr(xi   , n)
          lamp => dataptr(lam  , n)
          Ddp  => dataptr(Ddiag, n)

          lo = lwb(get_box(Q,n))
          hi = upb(get_box(Q,n))

          call orig_narrow_diffterm_3d(lo,hi,ng,dx,qp,upp,mup,xip,lamp,Ddp)

          call orig_hypterm_3d(lo,hi,ng,dx,up,qp,upp)
       end do

    end if

    wt2 = parallel_wtime()
    wt_hypdiff = wt_hypdiff + (wt2-wt1)

    if (update_courno) then
       call parallel_reduce(courno, courno_proc, MPI_MAX)
    end if

  end subroutine dUdt

  subroutine compute_courno(Q, dx, courno)
    type(multifab), intent(in) :: Q
    double precision, intent(in) :: dx(3)
    double precision, intent(inout) :: courno

    integer :: n, lo(3), hi(3), qlo(4), qhi(4)
    type(mfiter) :: mfi
    type(box) :: tbx
    double precision, pointer :: qp(:,:,:,:)

    !$omp parallel private(mfi, tbx, n, lo, hi, qlo, qhi, qp) &
    !$omp reduction(max:courno)

    call mfiter_build(mfi, Q, tiling=.true.)

    do while (next_tile(mfi,n))

       tbx = get_tilebox(mfi)
       lo = lwb(tbx)
       hi = upb(tbx)

       qp => dataptr(Q,n)
       qlo = lbound(qp)
       qhi = ubound(qp)

       call comp_courno_3d(lo,hi,dx,qp,qlo(1:3),qhi(1:3),courno)

    end do
    !$omp end parallel
  end subroutine compute_courno

end module advance_module


module adios2_adios

    use adios2_adios_init
    implicit none

contains

    subroutine adios2_declare_io(io, adios, io_name, ierr)
        integer(kind=8), intent(out) :: io
        integer(kind=8), intent(in) :: adios
        character*(*), intent(in)  :: io_name
        integer, intent(out) :: ierr

        call adios2_declare_io_f2c(io, adios, TRIM(ADJUSTL(io_name))//char(0), &
        & ierr)

    end

    subroutine adios2_finalize(adios, ierr)
        integer(kind=8), intent(in) :: adios
        integer, intent(out) :: ierr

        call adios2_finalize_f2c(adios, ierr)
    end

end module


       subroutine mgt_tillfactor(jj,bmix,emix,dtil,sol_thick)
	!!!!!!!!!!!!!!!!!!!!!!!
	! Armen 16 January 2008
	! This procedure increases tillage factor (tillagef(l,jj) per layer for each operation
	! The tillage factor settling will depend of soil moisture (tentatively) and must be called every day
	! For simplicity the settling is calculated now at the soil carbon subroutine because soil water content is available

	! The tillage factor depends on the cumulative soil disturbance rating = csdr
	! For simplicity, csdr is a function of emix
	! First step is to calculate "current" csdr by inverting tillage factor function
	! The effect of texture on tillage factor (ZZ) is removed first (and recovered at the end of the procedure)
	! YY = tillagef(l,jj) / ZZ
	! Since the tillage factor function is non linear, iterations are needed
	! XX = 0.5 is the initial value that works OK for the range of values observed
	! If a layer is only partially tilled then emix is corrected accordingly

	use soil_module

    implicit none

	integer, intent (in) :: jj        !none           |HRU number
    real, intent (in) :: bmix         !none           |biological mixing efficiency: this
                                      !               |number is zero for tillage operations
    integer :: l                      !none           |counter
    integer ::m1                      !none           |array location (see definition of ndays)
    integer :: m2                     !               |
    real :: emix                      !none           |mixing efficiency
    real :: dtil                      !mm             |depth of mixing
    real :: XX                        !varies         |variable to hold calculation results
	real :: sol_thick(soil(jj)%nly)   !               |
    integer :: j                      !none           |counter
    real :: zz                        !               |
    real :: yy                        !               |
    real :: xx1                       !               |
    real :: xx2                       !               |
    real :: csdr                      !               |

	emix = emix - bmix ! this is to avoid affecting tillage factor with biological mixing

	if (emix > 0.) then

	  do l = 1, soil(j)%nly

	    if (soil(jj)%phys(l)%d <= dtil) then
		  emix = emix
          else if(soil(jj)%phys(l)%d > dtil .AND. soil(jj)%phys(l-1)%d    &
                 < dtil) then
		emix = emix * (dtil - soil(jj)%phys(l-1)%d) / sol_thick(l)
	    else
		emix = 0.
	    end if

	    ! to save computation time if emix = 0 here then the other layers can be avoided
	    ! tillage always proceeds from top to bottom
	    if (emix == 0.) exit

	    xx = 0.
	    zz = 3. + (8. - 3.)*exp(-5.5*soil(jj)%phys(1)%clay/100.)
	    yy = soil(jj)%ly(l)%tillagef / zz
	    m1 = 1
	    m2 = 2

	    ! empirical solution for x when y is known and y=x/(x+exp(m1-m2*x))
	    if (yy > 0.01) then
		 xx1 = yy ** exp(-0.13 + 1.06 * yy)
		 xx2 = exp(0.64 + 0.64 * yy ** 100.)
		 xx = xx1 * xx2
	    end if

	    csdr = xx + emix
	    soil(jj)%ly(l)%tillagef = zz * (csdr / (csdr + exp(m1 - m2*csdr)))

	  end do

	end if

	return
	end subroutine mgt_tillfactor


text:
!===============================================================================
!> $Id: a9dba85759775b7901a5fe028b52124776dca40f $
!> Status bits for tasks
!===============================================================================
MODULE bits_mod
  integer, parameter::       static_bit=1
  integer, parameter::         leaf_bit=2
  integer, parameter::        ready_bit=2**2
  integer, parameter::     internal_bit=2**3
  integer, parameter::     boundary_bit=2**4
  integer, parameter::      virtual_bit=2**5
  integer, parameter::     external_bit=2**6
  integer, parameter:: swap_request_bit=2**7
  integer, parameter::      support_bit=2**8
  integer, parameter::    root_grid_bit=2**9
  integer, parameter::         star_bit=2**10
  integer, parameter::       planet_bit=2**11
  integer, parameter::         disk_bit=2**12
  integer, parameter::     particle_bit=2**13
  integer, parameter::        retry_bit=2**14
  integer, parameter::       remove_bit=2**15
  integer, parameter::   top_bottom_bit=2**16
! NOTE: the first (or last) 16 bits are masked out in some contexts; please CHECK
  integer, parameter::       has_BC_bit=2**30
  integer, parameter::        trace_bit=2**29
  integer, parameter::       frozen_bit=2**28
  integer, parameter::         busy_bit=2**27
  integer, parameter::  no_download_bit=2**26
  integer, parameter::   no_density_bit=2**25
  integer, parameter::        no_rt_bit=2**24
  integer, parameter::        nbors_bit=2**23   ! trigger init_nbors on vnbor
  integer, parameter::      counter_bit=2**20
  type bits_t
    integer:: remove       = remove_bit
    integer:: retry        = retry_bit
    integer:: ready        = ready_bit
    integer:: internal     = internal_bit
    integer:: virtual      = virtual_bit
    integer:: boundary     = boundary_bit
    integer:: external     = external_bit
    integer:: root_grid    = root_grid_bit
    integer:: not_leaf     = leaf_bit
    integer:: frozen       = frozen_bit
    integer:: swap_request = swap_request_bit
    integer:: support      = support_bit
    integer:: particle     = particle_bit
    integer:: star         = star_bit
    integer:: planet       = planet_bit
    integer:: static       = static_bit
    integer:: top_bottom   = top_bottom_bit
    integer:: has_BC       = has_BC_bit
    integer:: trace        = trace_bit
    integer:: busy         = busy_bit
    integer:: no_download  = no_download_bit
    integer:: no_density   = no_density_bit
    integer:: no_rt        = no_rt_bit
    integer:: init_nbors   = nbors_bit
    integer:: has_counter  = counter_bit
    integer:: hands_off    = ready_bit+busy_bit+virtual_bit+frozen_bit+external_bit
    !-------- BOUNDARY CONDITION BITS ------------------------------------------
    integer:: xl           = 2**1
    integer:: xu           = 2**2
    integer:: yl           = 2**3
    integer:: yu           = 2**4
    integer:: zl           = 2**5
    integer:: zu           = 2**6
    integer:: spherical    = 2**7
    integer:: symmetric    = 2**8
    integer:: antisymmetric= 2**9
    integer:: extrapolate  = 2**10
  end type
  type(bits_t):: bits
END MODULE bits_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
20091229 14:17:59 nbsp Hey, Id normally welcome you to the wiki, etc etc, but if you were surfing RC, it doesnt seem likely that youre new to the wiki.  Still, thanks for the recent edits: the wiki is built of contributions from everybody.  Since Woodland Cleaners makes the claim of being organic, I reinstated it, reworded it to say they are making a claim of, and then dropped in the rather important points you made on the talk page.  Since they are making the claim, it seems to be more informative and useful to the community to point out the issues with the claim than simply ignore it.  By the way  I saw your edit on the LARP entry; do you play? Users/JabberWokky Evan JabberWokky Edwards

  20091229 11:44:47   Howdy JW. Im new to daviswiki, but not wikis in general... Ive been checking this place out for a while, RC is a great place to get local news. Before I removed that businesses organic claim, I looked them up to see if theyre listed by CCOF as an organic business... no dice. My guess is they just threw it up as a marketing word and dont understand how involved it is getting organic certification, or the criminal and civil penalties involved in falsely claiming to be organic or certified organic. I used to play in LARPs, but back in their heyday in the late 90s and early oughts. It seems that WW hasnt been doing a good job of getting new players. Nice to meetcha. —EdwardNiemand

20091229 22:49:31 nbsp Yes, thanks for those comments about dry cleaners.  I found them very helpful. Users/CovertProfessor

20091230 16:43:01 nbsp Thanks for the tank house picture! Users/PhilipNeustrom

20100113 11:42:13 nbsp If youre doing a publicity campaign, want some Wiki Sticker/Requests stickers? Users/WilliamLewis
      Hey, those look dandy. I might take you up on the offer, but for now Im going to tentatively just flier... Id hate to get arrested for vandalism. Maybe if we start tabling again at the farmers market we could pass those out.  Users/EdwardNiemand

20100113 11:43:11 nbsp Ha!  Go for it... once somebody leads the way, you can likely get some support.  There have also been some other promotional things like tabling at the Farmers Market and other events. Users/JabberWokky

20100113 12:06:59 nbsp Hey, nice posters! They do a good job of playing up the wiki as a source of information. Do you have any good ideas to encourage contributions that arent reviews of businesses and such? Id like to get away from our perception as some sort of local, nonprofit Yelp (we have too many business owners interacting with us as such). Users/WilliamLewis

20100113 14:35:38 nbsp dont sign things as Ed <.<  It took me a few years to break that habit here :P Users/EdWins

20100113 20:54:44 nbsp Are you always a white knight?

Welcome to the wiki btw ♥ Users/StevenDaubert

20100115 12:51:03 nbsp @ Steven  Give me your weak, your downtrodden, your comic book shops with nasty comments, your dry cleaners yearning to breathe free... Users/EdwardNiemand

20100121 08:02:25 nbsp Im sorry its so long, there are several issues, and several board members I want to step down. Users/PxlAted

20100121 08:39:11 nbsp RE: I understand that youre pissed at them, and your position of homelessness sucks, but is bringing that feud onto this site what you should be spending your time on? You dont know me at all or what I spend my time on. Im not feuding, Im expressing specific problems, concerns, and complaints about certain members of a local, public, not for profit organization. Users/PxlAted

20100121 08:58:33 nbsp Personal problems and complaints. Users/EdwardNiemand

20100121 10:15:55 nbsp Wes and Dan are excluded from voting here because their own interests blind them to the greater good of the wiki. Do you have the power to do that? How do I get the power to prevent other people from voting? Can I declare your votes invalid, too?  As a side comment: Well said.  I actually dropped the phrase, but I originally wrote too much pseudolegalism, authoritarian posturing and game playing.  I figured that would obscure the more ontopic points, and many editors present their opinions with great assertiveness and might have thought I was referring to them.  I think William was just trying to move things along and block potential complaints, aware that anybody who wanted to could object. Users/JabberWokky

   Just making sure that this place isnt succumbing to wikidisease. EN

20100128 10:40:46 nbsp Our IRS tax exemption prohibits pricelists from forprofit entities. When I removed the prices on the Real Computers page, Brian fought with me, found out I was right, and started removing prices from all his competitors. Users/WilliamLewis

20100407 16:12:54 nbsp http://daviswiki.org/Con_Artists?actiondiff&version2161&version1160 I was going to say, cpc § 148.5(b). Users/WilliamLewis

20100407 16:14:06 nbsp For someone who was worried about the Con Artist page being too cluttered, you sure have gone off on a tangent. Users/CovertProfessor

      Ha ha! I know, I know, I knew it when I was writing it, but its hard to be both right and quiet about it! :) I moved the thread down into the comments section of the page, to unclutter it a bit. Users/EdwardNiemand

20100407 16:52:43 nbsp More important than pricing information is the phrasing of the page.  That page looks like its written in the tone of the business we should rewrite it Users/PhilipNeustrom
    Agreed. EN

20100623 21:08:22 nbsp Can you send me an email?  philipn at gmail.com Users/PhilipNeustrom

20100624 09:29:48 nbsp email sent. Users/EdwardNiemand

20100627 16:56:18 nbsp Yup, thats me in the video.   It got me to install the flash plugin that I had been stubbornly not plugging in, just to see my eyes widening into interrobangs. Users/Angel.York

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE GSGMGN  ( xl, yb, xr, yt, iret )
C************************************************************************
C* GSGMGN								*
C* 									*
C* This subroutine sets the margin sizes to be used in graph mode of 	*
C* the map/graph coordinate system.  The margin sizes may be specified	*
C* as either a fraction of the view region or as a multiple of the	*
C* current character size.  If the values entered are greater than 0    *
C* and less than 1, they are considered to be a fraction of the view    *
C* region.  If the values are 1 or greater, they are taken to be        *
C* multiples of the character size at the time of the call to GSGMGN.   *
C* The default margin size is zero.					*
C* 									*
C* GSGMGN should be called before any plotting is done; margin size	*
C* should not be changed after plotting has begun.  If size is 		*
C* specified in terms of character size, the text size at the time 	*
C* of the call is used.  Later changes to the size do not affect 	*
C* the margins.								*
C* 									*
C* GSGMGN  ( XL, YB, XR, YT, IRET )					*
C*									*
C* Input parameters:							*
C* 	XL		REAL		Left margin size		*
C* 	YB		REAL		Bottom margin size		*
C* 	XR		REAL		Right margin size		*
C* 	YT		REAL		Top margin size 		*
C* 									*
C* Output parameters:							*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* I. Graffman/RDS	 5/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 9/88	GEMPAK4					*
C* L. Williams/EAi       3/94   Removed blank comments from header      *
C* A. Hardy/GSC          6/98   Cleaned up prolog                       *
C************************************************************************
	INCLUDE 	'FUNCCODE.PRM'
	INCLUDE 	'ERROR.PRM'
C*
	INTEGER 	isend (2)
	REAL 		rsend (4)
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 6
	isend (2) = FSGMGN
	rsend (1) = xl
	rsend (2) = yb
	rsend (3) = xr
	rsend (4) = yt
	CALL GPUT  ( isend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
	CALL GPUTR  ( rsend, 4, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
C*	Get output parameters.
C
	CALL GGET  ( iret, 1, ier )
	IF  ( ier .ne. NORMAL )  iret = ier
C*
	RETURN
	END


text:
      SUBROUTINE PARPT (RCX,RCY, PARITY, PRVC1X,PRVC1Y, PRVC2X,PRVC2Y,
     1                  XA,YA, X,Y, STEP1,STEP2)

C  PARPT accepts a reference point and 2 vectors and computes a
C  corner of the smallest parallelogram containing a specified
C  rectangle. The corner is found in the direction specified by
C  parity.

C  Input parameters
C    RCX,RCY       : reference point within rectangle
C    PARITY        : +1 to search in direction of vectors, -1 to
C                    search in opposite direction
C    PRVC1X,PRVC1Y : first vector to increment reference point
C    PRVC2X,PRVC2Y : second vector to increment reference point
C    XA,YA         : 4 corners of rectangle to be contained
C  Output parameters
C    X,Y           : corner of parallelogram
C    STEP1,STEP2   : number of steps taken from reference point in
C                    direction of 1st,2nd vector

      INTEGER  PARITY, STEP1, STEP2

      REAL     RCX,RCY, PRVC1X,PRVC1Y, PRVC2X,PRVC2Y, XA(4),YA(4), X,Y
      REAL     A1,B1,C1, A2,B2,C2

C  find first, second boundary lines
      CALL PARLIN (RCX,RCY, PARITY, PRVC1X,PRVC1Y, PRVC2X,PRVC2Y,
     1             XA,YA, A1,B1,C1, STEP1)
      CALL PARLIN (RCX,RCY, PARITY, PRVC2X,PRVC2Y, PRVC1X,PRVC1Y,
     1             XA,YA, A2,B2,C2, STEP2)
C  intersection of lines is corner
      CALL LINTPT (A1,B1,C1, A2,B2,C2, X,Y)

      END


text:
module mod_fluid
#include <messenger.h>
    use mod_kinds,              only: rk, ik
    use mod_constants,          only: ZERO, ONE
    implicit none

    ! Fluid equation constants
    real(rk), parameter :: Rgas  = 287.15_rk            ! Specific gas constant
    real(rk), parameter :: gam   = 1.4_rk               ! Ratio of specific heats
    real(rk), parameter :: cp    = Rgas*gam/(gam-ONE)   ! Specific heat capacity, constant pressure
    real(rk)            :: omega = 0._rk                ! Rotation rate (rad/s)
    real(rk)            :: Pr    = 0.72_rk              ! Prandtl number

end module mod_fluid

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
      SUBROUTINE FGAUSS(X,A,Y,DYDA,NA)
      DIMENSION A(NA),DYDA(NA)
      Y=0.
      DO 11 I=1,NA-1,3
        ARG=(X-A(I+1))/A(I+2)
        EX=EXP(-ARG**2)
        FAC=A(I)*EX*2.*ARG
        Y=Y+A(I)*EX
        DYDA(I)=EX
        DYDA(I+1)=FAC/A(I+2)
        DYDA(I+2)=FAC*ARG/A(I+2)
11    CONTINUE
      RETURN
      END

ext: for


text:
C    %W% %G%
      subroutine govint(ign, edt, iabort)
c     this routine initializes the governors

      implicit none

      include 'tspinc/params.inc'
      include 'tspinc/machd1.inc'
      include 'tspinc/gov.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/spare1.inc'

c     local variables
c     ign    -  generator number. i and j are local counting variables.
c     mgov   -  governor type, igv is the govdat offset
c     gvpwr  -  power out of governor (govpwr contains PM1 for cc govs)
c     iabort -  flag to abort run
c     x1, x2, x3p, x3 -> x14       states

C     Local variables

      integer ign, i, j, mgov, igv, iabort
      real    edt, dt, x1, x2, x3p, x3, x4, x5, x6, x7, x8,
     &        x9, x10, x11, x12, x13, x14, gvpwr, xtemp

      dt = edt/60.0

c     assign offset to govdat array for local use
      igv = govindx(ign)

      mgov = govtyp(ign)

C     assign large value to govmx
      govmx(ign) = 1.0e6

c     zero govpwr if very small
      if (abs(govpwr(ign)) .lt. 1.0e-6) govpwr(ign) = 0.0

c     determine position of switch in feedback path and save
      if (govdat(igv+16).lt.1.0e-6) then
        govmode(ign) = 1
      else
        govmode(ign) = 2
      endif

c     initialize turbine
      call turbint(ign, govpwr(ign), gvpwr, dt)

c     initialize nonlinear curve
      if (govdat(igv+15).ne.0.0.and.govtyp(ign).ne.5) then
        call cvint( govdat(igv+15), govcv(ign), gvpwr, xtemp,
     &              govdat(igv+ 1) )
        gvpwr = xtemp
      endif

c     put govpwr into Phigh if Phigh is zero and into Plow if Plow is zero
      if (govdat(igv+ 9).lt.1.0e-6) govdat(igv+9) = gvpwr
      if (govdat(igv+10).lt.1.0e-6) govdat(igv+10) = gvpwr

      if (mgov.eq.1) then
c ********************* GS governor *******************************

c       store initial governor power
        govdat(igv+17) = gvpwr

c       store time constants in convenient format
        govdat(igv+ 4) = govdat(igv+ 4)*2.0/dt + 1.0
        govdat(igv+ 5) = govdat(igv+ 5)*2.0/dt + 1.0
        govdat(igv+ 6) = 1.0/govdat(igv+ 6)
        govdat(igv+16) = govdat(igv+16)*2.0/dt + 1.0

c       initialize states
        x1  = 0.0
        x2  = 0.0
        x3p = 0.0
        x3  = govdat(igv+17)
        x4  = x3
        x5  = x4
        x6  = x5

c       store states that need to be stored
        govdat(igv+21) = x4
        govdat(igv+22) = x5
        govdat(igv+24) = 0.0
        govdat(igv+25) = x1
        govdat(igv+26) = x3
        govdat(igv+27) = x4

c       assign history values
        govdat(igv+18) =   govdat(igv+3)*(govdat(igv+4) - 2.0)*x2
     &                                  -(govdat(igv+5) - 2.0)*x1
        govdat(igv+19) = x3*2.0/dt + x3p
        govdat(igv+20) = x6*(govdat(igv+16)-2.0) + gvpwr

      else if (mgov.eq.2) then
c ********************* GH governor *******************************

c       store initial governor power * steady state droop
        govdat(igv+17) = gvpwr*govdat(igv+3)

c       store time constants in convenient format
        govdat(igv+ 5) = govdat(igv+ 5)*2.0/dt + 1.0
        govdat(igv+19) = 2.0*govdat(igv+ 6)*govdat(igv+18)/dt
        govdat(igv+ 6) = govdat(igv+ 6)*2.0/dt + 1.0
        govdat(igv+16) = govdat(igv+16)*2.0/dt + 1.0

c       initialize states
        x1  = 0.0
        x2  = 0.0
        x3  = 0.0
        x4  = gvpwr
        x5  = x4
        x6  = x5
        x7  = 0.0
        x8  = x5

c       store states that need to be stored
        govdat(igv+24) = x5
        govdat(igv+25) = x8
        govdat(igv+30) = x1
        govdat(igv+31) = x4
        govdat(igv+32) = x5

c       assign history values
        govdat(igv+20) = govdat(igv+4)*(govdat(igv+5)-1.0)*x3 + x2
        govdat(igv+21) = x4*2.0/dt + x3
        govdat(igv+22) = (govdat(igv+6)-1.0)*x7 - govdat(igv+19)*x5
        govdat(igv+23) = x6*(govdat(igv+16)-2.0) + gvpwr

      else if (mgov.eq.3) then
c ********************* GP governor *******************************

c       store initial governor power * steady state droop
        govdat(igv+17) = gvpwr*govdat(igv+3)

c       store time constants in convenient format
        govdat(igv+ 4) = govdat(igv+ 4)*2./dt + 1.0
        govdat(igv+ 5) = govdat(igv+ 5)*2./dt + 1.0
        govdat(igv+ 6) = govdat(igv+ 6)*2./dt + 1.0
        govdat(igv+16) = govdat(igv+16)*2./dt + 1.0

c       initialize states
        x1  = 0.0
        x2  = 0.0
        x3  = 0.0
        x4  = gvpwr
        x5  = 0.0
        x6  = gvpwr
        x7  = gvpwr
        x8  = 0.0
        x9  = 0.0
        x10 = gvpwr
        x11 = gvpwr
        x12 = gvpwr

c       store states that need to be stored
        govdat(igv+32) = x6
        govdat(igv+33) = x11
        govdat(igv+34) = x12
        govdat(igv+37) = x1
        govdat(igv+38) = x10
        govdat(igv+39) = x11

c       assign history values
        govdat(igv+22) = x2 + x3*(govdat(igv+4)-2.0)
        govdat(igv+23) = x3*govdat(igv+20) + x4*2.0/dt
        govdat(igv+24) = x5*(govdat(igv+5)-2.) - x3*2.*govdat(igv+19)/dt
        govdat(igv+25) = x8*govdat(igv+21) + x9*(govdat(igv+6)-2.)
        govdat(igv+27) = x9 + x10*2.0/dt
        govdat(igv+35) = x7*(govdat(igv+16)-2.0) + gvpwr

      else if (mgov.eq.4) then
c ********************* G2 governor *******************************

c       store initial governor power * steady state droop
        govdat(igv+17) = gvpwr*govdat(igv+3)

c       store time constants in convenient format
        govdat(igv+ 4) = govdat(igv+ 4)*2./dt + 1.0
        govdat(igv+ 5) = govdat(igv+ 5)*2./dt + 1.0
        govdat(igv+ 6) = govdat(igv+ 6)*2./dt + 1.0
        govdat(igv+16) = govdat(igv+16)*2./dt + 1.0

c       initialize states
        x1  = 0.0
        x2  = 0.0
        x3  = 0.0
        x4  = 0.0
        x5  = 0.0
        x6  = 0.0
        x7  = 0.0
        x8  = gvpwr
        x9  = gvpwr
        x10 = 0.0
        x11 = 0.0
        x12 = gvpwr
        x13 = gvpwr
        x14 = gvpwr

c       store states that need to be stored
        govdat(igv+34) = x6
        govdat(igv+35) = x8
        govdat(igv+36) = x13
        govdat(igv+37) = x14
        govdat(igv+40) = x1
        govdat(igv+41) = x12
        govdat(igv+42) = x13

c       assign history values
        govdat(igv+22) = x2*(govdat(igv+4)-2.0) + x1
        govdat(igv+23) = x3*(govdat(igv+5)-2.0)
     &                   - 2.0/dt*govdat(igv+18)*x2
        govdat(igv+24) = x4*(govdat(igv+5)-2.0)
     &                   - 2.0/dt*govdat(igv+19)*x2
        govdat(igv+25) = x5*(govdat(igv+5)-2.0)
     &                   - 2.0/dt*x4
        govdat(igv+27) = 2.0/dt*x8 + govdat(igv+20)*x7
        govdat(igv+28) = x11*(govdat(igv+6)-2.0) + govdat(igv+21)*x10
        govdat(igv+29) = 2.0/dt*x12 + x11
        govdat(igv+38) = x9*(govdat(igv+16)-2.0) + govdat(igv+26)

      else if (mgov.eq.5) then
c ********************* GW governor *******************************

c       store initial governor power
        govdat(igv+7) = gvpwr

c       store time constants in convenient format
        govdat(igv+4) = 2.0*govdat(igv+4)/dt + 1.0
        govdat(igv+5) = 2.0*govdat(igv+5)/dt + 1.0
        govdat(igv+6) = 2.0*govdat(igv+6)/dt + 1.0

c       initialize states
        x1 = 0.0
        x2 = 0.0
        x3 = 0.0
        x4 = gvpwr

c       assign history values
        govdat(igv+8) = x2*(govdat(igv+5)-2.0) - x1*(govdat(igv+4)-2.0)
        govdat(igv+9) = x3*(govdat(igv+6)-2.0) + x2

      endif

c       test for governor power exceeding limits
        if (gvpwr .gt. 1.001*govdat(igv+1) .or.
     &      gvpwr .lt. 0.999*govdat(igv+2)) then
          write(errbuf(1),100) govname(ign), govdat(igv+2),
     &                         gvpwr, govdat(igv+1)
100       format('Initial power at ', a8, ' exceeds governor limits: ',
     &           f7.3, ' <', f7.3, ' <', f7.3)
          call prterr('E', 1)
          iabort = 1
        endif

c     Put Phigh into Pmax and/or Plow into Pmin depending on Iblock
      if (abs(govdat(igv+11)-1.0).lt.1.0e-6) then
        govdat(igv+2) = govdat(igv+10)
      else if (abs(govdat(igv+11)-2.0).lt.1.0e-6) then
        govdat(igv+1) = govdat(igv+ 9)
      else if (abs(govdat(igv+11)-3.0).lt.1.0e-6) then
        govdat(igv+2) = govdat(igv+10)
        govdat(igv+1) = govdat(igv+ 9)
      endif

      return
      end


text:
The Davis Manor neighborhood is located between 5th Street and 8th Street, and between L Street and Pole Line Road (as defined on their http://members.dcn.org/dmnc/index.html website). The Davis Manor Center Shopping Centers shopping center is located in this neighborhood and its revival an ongoing issue for the Neighborhood Council.

The Manor Neighborhood is basically East Davis, or old east Davisbefore Mace Ranch. The green belt the encircles Mace Ranch can be seen as the old limits of east Davis. There is the Manor shopping center and Manor park, which has a public swimming pool and a huge cement slide. I am not sure what the extent of the Manor Neighborhood is. This maybe a more instituted form of what is also know as the Peoples Republic of Davis Peoples Republic of East Davis (PRED).

Davis Manor encompasses much more than that.  The area surrounding Slide Hill Park (Manor Park) also is part of the neighborhood.  Users/TobinJones

20080521 22:56:03 nbsp Davis Manor was originally the name of a housing development, and the DMNC is for that area. I respect that Tobin has a different view, but I think that since this page is about the Council, the content about Manor Pool, West Manor, and the Birch Lane/Tulip Lane area is offtopic enough to perhaps warrant future deletion. Users/DougWalter