text:
C
C $Id: ngdots.f,v 1.14 2008-07-27 00:17:17 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE NGDOTS(X,Y,NUM,SIZE,ICOLOR)
C
C  Draws filled dots at positions (X(I),Y(I),I=1,NUM) at size
C  SIZE with color given by the color index ICOLOR.  Size is
C  the diameter of the dot in world coordinates along the vertical
C  axis.
C
C  The dots are scaled appropriately so that they will be circular
C  when the normalization transformation does not preserve aspect
C  ratio.
C
C  The algorithm is to construct a single circle and then translate
C  it to the various coordinate positions.  The original circle is
C  computed by using trig functions to get points in the first
C  octant, and then using symmetries to get the rest of the points
C  on the circle.
C
C  The number of points used for the circle is adjusted depending
C  on the relative size of the circle.  The maximum number of points
C  is 512 and the minimum is 8.
C
      SAVE
C
      include 'ngcomn.h'
C
      DIMENSION X(NUM),Y(NUM),TWIN(4),TVPT(4)
      DIMENSION CIRCX(513),CIRCY(513),CIRCXX(513),CIRCYY(513)
      INTEGER   OLDASF(13),NEWASF(13)
C
C  PI and number of points per unit of arc length.
C
      DATA PI,NPU/3.1415927,1200/
C
C  Size below which markers are drawn in addition to the filled dots
C  to insure that the dots do not disappear by falling between pixels.
C
      DATA DM/.05/
C
C  Size cutoff for increasing the number of points at smaller sizes
C  so that the dots look round.
C
      DATA PC/.03/
C
C  Minimum dot size.  Below this tolerance, only dot markers are drawn.
C
      DATA SM/.001/
C
      DATA IFIRST/0/
C
C  If log scaling or mirror imaging has been requested, issue a
C  warning that only GKS world coordinates are accepted.
C
      CALL GETUSV('LS',LSV)
      CALL GETUSV('MI',MIV)
      IF ((LSV.NE.1 .OR. MIV.NE.1) .AND. IFIRST.EQ.0) THEN
        PRINT *, 'NGDOTS -- **Warning**'
        PRINT *, '          Log scaling or mirror-imaging have been requ
     +ested via a SET call;'
        PRINT *, '          NGDOTS accepts only GKS world coordinates.'
        IFIRST = 1
      ENDIF
C
C  Calculate the number of points to use in calculating the points
C  in the first octant of the circle.
C
      CALL GQCNTN(IER,NTNR)
      CALL GQNT(NTNR,IER,TWIN,TVPT)
C
C  Normalize the size.
C
      SD = SIZE*(TVPT(4)-TVPT(3))/(TWIN(4)-TWIN(3))
C
C  Set a minimum dot size.
C
      IF (SD .LT. SM) THEN
        SD = SM
        TSIZE = SD*(TWIN(4)-TWIN(3))/(TVPT(4)-TVPT(3))
      ELSE
        TSIZE = SIZE
      ENDIF
C
C  Calculate the number of points on the curve.
C
      NP = INT(NPU*PI*SD)
C
C  Increase the number of points for smaller circles.
C
      IF (SD .LE. PC) NP = NP*2
      NP = NP/8
      NP = NP+1
      NP = MAX(NP,2)
      NP = MIN(NP,64)
C
C  Center the dot at (0.,0.) and translate it to the (X,Y)
C  coordinates later.
C
C  Calculate the points in the first octant.
C
      AINC = .25*PI/REAL(NP)
      R = .5*TSIZE
      DO 10 I=1,NP
        T = REAL(I-1)*AINC
        CIRCX(I) = R*COS(T)
        CIRCY(I) = R*SIN(T)
   10 CONTINUE
C
C  Use symmetries to get the other points on the circle.
C
      NP2 = 2*NP
      DO 30 I=NP+1,NP2-1
        CIRCX(I) = CIRCY(NP2-I)
        CIRCY(I) = CIRCX(NP2-I)
   30 CONTINUE
      NP = NP2-1
      NP2 = 2*NP
      DO 40 I=NP+1,NP2-1
        CIRCX(I) = -CIRCX(NP2-I)
        CIRCY(I) =  CIRCY(NP2-I)
   40 CONTINUE
      NP = NP2-1
      NP2 = 2*NP
      DO 50 I=NP+1,NP2-1
        CIRCX(I) =  CIRCX(NP2-I)
        CIRCY(I) = -CIRCY(NP2-I)
   50 CONTINUE
      NP = NP2-1
C
C  Scale factor for non-square aspect ratios.
C
      WSCALE = (TWIN(2)-TWIN(1))/(TWIN(4)-TWIN(3))
      VSCALE = (TVPT(4)-TVPT(3))/(TVPT(2)-TVPT(1))
      DO 20 I=1,NP
        CIRCX(I) = WSCALE*VSCALE*CIRCX(I)
   20 CONTINUE
C
C  Draw filled dots if the size is large enough.
C
      IF (SD .GE. SM) THEN
C
C  Save the fill color index and interior style.
C
        CALL GQFACI(IER,ICOLRO)
        CALL GQFAIS(IER,ISTYLO)
C
C  Set the fill style to solid and the fill collor to ICOLOR.
C
        CALL GSFAIS(1)
        CALL GSFACI(ICOLOR)
C
C  Draw the dots or circles, depending on ICDFLG.
C
        DO 60 I=1,NUM
          DO 70 J=1,NP
            CIRCXX(J) = X(I) + CIRCX(J)
            CIRCYY(J) = Y(I) + CIRCY(J)
   70     CONTINUE
          IF (ICDFLG .EQ. 0) THEN
            CALL GFA(NP,CIRCXX,CIRCYY)
          ELSE
            CALL GQPLCI(IER,IPLCO)
            CALL GSPLCI(ICOLOR)
            CALL GPL(NP,CIRCXX,CIRCYY)
            CALL GSPLCI(IPLCO)
          ENDIF
   60   CONTINUE
C
C  Restore the fill color and style.
C
        CALL GSFACI(ICOLRO)
        CALL GSFAIS(ISTYLO)
      ENDIF
C
C  Draw polymarkers if the size is small enough (note that there is
C  an overlap where both dots and markers are drawn).  The polymarkers
C  are drawn only when filled dots are requested and not when circles
C  are requested.
C
      IF (SD.LT.DM .AND. ICDFLG.EQ.0) THEN
C
C  Save current GKS attributes.
C
        CALL GQASF(IER,OLDASF)
        CALL GQMK (IER,MTYPEO)
        CALL GQMKSC(IER,SZSFO)
        CALL GQPMCI(IER,IMCOLO)
C
C  Set attributes for polymarkers and draw them.
C
        DO 80 I=1,13
          NEWASF(I) = OLDASF(I)
   80   CONTINUE
        NEWASF(4) = 1
        NEWASF(5) = 1
        NEWASF(6) = 1
        CALL GSASF(NEWASF)
        CALL GSMK(1)
        CALL GSMKSC(1.0)
        CALL GSPMCI(ICOLOR)
        CALL GPM(NUM,X,Y)
C
C  Restore original attributes.
C
        CALL GSASF(OLDASF)
        CALL GSMK(MTYPEO)
        CALL GSMKSC(SZSFO)
        CALL GSPMCI(IMCOLO)
      ENDIF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
C Modified by J. Vandenplas (2019/10/28) to allow semi-positive definite
C matrix
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
C
C  =========== DOCUMENTATION ===========
C
C Online html documentation available at
C            http://www.netlib.org/lapack/explore-html/
C
C  Definition:
C  ===========
C
C       SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
C
C       .. Scalar Arguments ..
C       DOUBLE PRECISION ALPHA
C       INTEGER LDA,LDB,M,N
C       CHARACTER DIAG,SIDE,TRANSA,UPLO
C       ..
C       .. Array Arguments ..
C       DOUBLE PRECISION A(LDA,*),B(LDB,*)
C       ..
C
C
C> \par Purpose:
C  =============
C>
C> \verbatim
C>
C> DTRSM  solves one of the matrix equations
C>
C>    op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
C>
C> where alpha is a scalar, X and B are m by n matrices, A is a unit,
Cor
C> non-unit,  upper or lower triangular matrix  and  op( A )  is one
Cof
C>
C>    op( A ) = A   or   op( A ) = A**T.
C>
C> The matrix X is overwritten on B.
C> \endverbatim
C
C  Arguments:
C  ==========
C
C> \param[in] SIDE
C> \verbatim
C>          SIDE is CHARACTER*1
C>           On entry, SIDE specifies whether op( A ) appears on the
Cleft
C>           or right of X as follows:
C>
C>              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
C>
C>              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
C> \endverbatim
C>
C> \param[in] UPLO
C> \verbatim
C>          UPLO is CHARACTER*1
C>           On entry, UPLO specifies whether the matrix A is an upper
Cor
C>           lower triangular matrix as follows:
C>
C>              UPLO = 'U' or 'u'   A is an upper triangular matrix.
C>
C>              UPLO = 'L' or 'l'   A is a lower triangular matrix.
C> \endverbatim
C>
C> \param[in] TRANSA
C> \verbatim
C>          TRANSA is CHARACTER*1
C>           On entry, TRANSA specifies the form of op( A ) to be used
Cin
C>           the matrix multiplication as follows:
C>
C>              TRANSA = 'N' or 'n'   op( A ) = A.
C>
C>              TRANSA = 'T' or 't'   op( A ) = A**T.
C>
C>              TRANSA = 'C' or 'c'   op( A ) = A**T.
C> \endverbatim
C>
C> \param[in] DIAG
C> \verbatim
C>          DIAG is CHARACTER*1
C>           On entry, DIAG specifies whether or not A is unit
Ctriangular
C>           as follows:
C>
C>              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
C>
C>              DIAG = 'N' or 'n'   A is not assumed to be unit
C>                                  triangular.
C> \endverbatim
C>
C> \param[in] M
C> \verbatim
C>          M is INTEGER
C>           On entry, M specifies the number of rows of B. M must be
Cat
C>           least zero.
C> \endverbatim
C>
C> \param[in] N
C> \verbatim
C>          N is INTEGER
C>           On entry, N specifies the number of columns of B.  N must
Cbe
C>           at least zero.
C> \endverbatim
C>
C> \param[in] ALPHA
C> \verbatim
C>          ALPHA is DOUBLE PRECISION.
C>           On entry,  ALPHA specifies the scalar  alpha. When  alpha
Cis
C>           zero then  A is not referenced and  B need not be set
Cbefore
C>           entry.
C> \endverbatim
C>
C> \param[in] A
C> \verbatim
C>          A is DOUBLE PRECISION array, dimension ( LDA, k ),
C>           where k is m when SIDE = 'L' or 'l'
C>             and k is n when SIDE = 'R' or 'r'.
C>           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by
Ck
C>           upper triangular part of the array  A must contain the
Cupper
C>           triangular matrix  and the strictly lower triangular part
Cof
C>           A is not referenced.
C>           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by
Ck
C>           lower triangular part of the array  A must contain the
Clower
C>           triangular matrix  and the strictly upper triangular part
Cof
C>           A is not referenced.
C>           Note that when  DIAG = 'U' or 'u',  the diagonal elements
Cof
C>           A  are not referenced either,  but are assumed to be
Cunity.
C> \endverbatim
C>
C> \param[in] LDA
C> \verbatim
C>          LDA is INTEGER
C>           On entry, LDA specifies the first dimension of A as
Cdeclared
C>           in the calling (sub) program.  When  SIDE = 'L' or 'l'
Cthen
C>           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or
C'r'
C>           then LDA must be at least max( 1, n ).
C> \endverbatim
C>
C> \param[in,out] B
C> \verbatim
C>          B is DOUBLE PRECISION array, dimension ( LDB, N )
C>           Before entry,  the leading  m by n part of the array  B
Cmust
C>           contain  the  right-hand  side  matrix  B,  and  on exit
Cis
C>           overwritten by the solution matrix  X.
C> \endverbatim
C>
C> \param[in] LDB
C> \verbatim
C>          LDB is INTEGER
C>           On entry, LDB specifies the first dimension of B as
Cdeclared
C>           in  the  calling  (sub)  program.   LDB  must  be  at
Cleast
C>           max( 1, m ).
C> \endverbatim
C
C  Authors:
C  ========
C
C> \author Univ. of Tennessee
C> \author Univ. of California Berkeley
C> \author Univ. of Colorado Denver
C> \author NAG Ltd.
C
C> \date December 2016
C
C> \ingroup double_blas_level3
C
C> \par Further Details:
C  =====================
C>
C> \verbatim
C>
C>  Level 3 Blas routine.
C>
C>
C>  -- Written on 8-February-1989.
C>     Jack Dongarra, Argonne National Laboratory.
C>     Iain Duff, AERE Harwell.
C>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
C>     Sven Hammarling, Numerical Algorithms Group Ltd.
C> \endverbatim
C>
C  =====================================================================
      SUBROUTINE dgtrsm(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
C
C  -- Reference BLAS level3 routine (version 3.7.0) --
C  -- Reference BLAS is a software package provided by Univ. of
C  Tennessee,    --
C  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG
C  Ltd..--
C     December 2016
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER LDA,LDB,M,N
      CHARACTER DIAG,SIDE,TRANSA,UPLO
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(lda,*),B(ldb,*)
C     ..
C
C  =====================================================================
C
C     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
C     ..
C     .. External Subroutines ..
      EXTERNAL xerbla
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC max
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION atmp
      DOUBLE PRECISION TEMP
      INTEGER I,INFO,J,K,NROWA
      LOGICAL LSIDE,NOUNIT,UPPER
C     ..
C     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      DOUBLE PRECISION TOL
      parameter(one=1.0d+0,zero=0.0d+0)
      parameter(tol=1.0d-10)
C     ..
C
C     Test the input parameters.
C
      lside = lsame(side,'L')
      IF (lside) THEN
          nrowa = m
      ELSE
          nrowa = n
      END IF
      nounit = lsame(diag,'N')
      upper = lsame(uplo,'U')
C
      info = 0
      IF ((.NOT.lside) .AND. (.NOT.lsame(side,'R'))) THEN
          info = 1
      ELSE IF ((.NOT.upper) .AND. (.NOT.lsame(uplo,'L'))) THEN
          info = 2
      ELSE IF ((.NOT.lsame(transa,'N')) .AND.
     +         (.NOT.lsame(transa,'T')) .AND.
     +         (.NOT.lsame(transa,'C'))) THEN
          info = 3
      ELSE IF ((.NOT.lsame(diag,'U')) .AND. (.NOT.lsame(diag,'N')))THEN
          info = 4
      ELSE IF (m.LT.0) THEN
          info = 5
      ELSE IF (n.LT.0) THEN
          info = 6
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 9
      ELSE IF (ldb.LT.max(1,m)) THEN
          info = 11
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('DTRSM ',info)
          RETURN
      END IF
C
C     Quick return if possible.
C
      IF (m.EQ.0 .OR. n.EQ.0) RETURN
C
C     And when  alpha.eq.zero.
C
      IF (alpha.EQ.zero) THEN
          DO 20 j = 1,n
              DO 10 i = 1,m
                  b(i,j) = zero
   10         CONTINUE
   20     CONTINUE
          RETURN
      END IF
C
C     Start the operations.
C
      IF (lside) THEN
          IF (lsame(transa,'N')) THEN
C
C           Form  B := alpha*inv( A )*B.
C
              IF (upper) THEN
                  DO 60 j = 1,n
                      IF (alpha.NE.one) THEN
                          DO 30 i = 1,m
                              b(i,j) = alpha*b(i,j)
   30                     CONTINUE
                      END IF
                      DO 50 k = m,1,-1
                          IF (b(k,j).NE.zero) THEN
!                              IF (nounit) b(k,j) = b(k,j)/a(k,k)
                              IF (nounit)then
                               atmp=0.d0
                               if(a(k,k).gt.tol)atmp=1.d0/a(k,k)
                               b(k,j) = b(k,j)*atmp
                              endif
                              DO 40 i = 1,k - 1
                                  b(i,j) = b(i,j) - b(k,j)*a(i,k)
   40                         CONTINUE
                          END IF
   50                 CONTINUE
   60             CONTINUE
              ELSE
                  DO 100 j = 1,n
                      IF (alpha.NE.one) THEN
                          DO 70 i = 1,m
                              b(i,j) = alpha*b(i,j)
   70                     CONTINUE
                      END IF
                      DO 90 k = 1,m
                          IF (b(k,j).NE.zero) THEN
!                              IF (nounit) b(k,j) = b(k,j)/a(k,k)
                              IF (nounit)then
                               atmp=0.d0
                               if(a(k,k).gt.tol)atmp=1.d0/a(k,k)
                               b(k,j) = b(k,j)*atmp
                              endif
                              DO 80 i = k + 1,m
                                  b(i,j) = b(i,j) - b(k,j)*a(i,k)
   80                         CONTINUE
                          END IF
   90                 CONTINUE
  100             CONTINUE
              END IF
          ELSE
C
C           Form  B := alpha*inv( A**T )*B.
C
              IF (upper) THEN
                  DO 130 j = 1,n
                      DO 120 i = 1,m
                          temp = alpha*b(i,j)
                          DO 110 k = 1,i - 1
                              temp = temp - a(k,i)*b(k,j)
  110                     CONTINUE
!                          IF (nounit) temp = temp/a(i,i)
                          IF (nounit)then
                           atmp=0.d0
                           if(a(i,i).gt.tol)atmp=1.d0/a(i,i)
                           temp = temp*atmp
                          endif
                          b(i,j) = temp
  120                 CONTINUE
  130             CONTINUE
              ELSE
                  DO 160 j = 1,n
                      DO 150 i = m,1,-1
                          temp = alpha*b(i,j)
                          DO 140 k = i + 1,m
                              temp = temp - a(k,i)*b(k,j)
  140                     CONTINUE
!                          IF (nounit) temp = temp/a(i,i)
                          IF (nounit)then
                           atmp=0.d0
                           if(a(i,i).gt.tol)atmp=1.d0/a(i,i)
                           temp = temp*atmp
                          endif
                          b(i,j) = temp
  150                 CONTINUE
  160             CONTINUE
              END IF
          END IF
      ELSE
          IF (lsame(transa,'N')) THEN
C
C           Form  B := alpha*B*inv( A ).
C
              IF (upper) THEN
                  DO 210 j = 1,n
                      IF (alpha.NE.one) THEN
                          DO 170 i = 1,m
                              b(i,j) = alpha*b(i,j)
  170                     CONTINUE
                      END IF
                      DO 190 k = 1,j - 1
                          IF (a(k,j).NE.zero) THEN
                              DO 180 i = 1,m
                                  b(i,j) = b(i,j) - a(k,j)*b(i,k)
  180                         CONTINUE
                          END IF
  190                 CONTINUE
                      IF (nounit) THEN
!                          temp = one/a(j,j)
                          atmp=0.d0
                          if(a(j,j).gt.tol)atmp=1.d0/a(j,j)
                          temp = one*atmp
                          DO 200 i = 1,m
                              b(i,j) = temp*b(i,j)
  200                     CONTINUE
                      END IF
  210             CONTINUE
              ELSE
                  DO 260 j = n,1,-1
                      IF (alpha.NE.one) THEN
                          DO 220 i = 1,m
                              b(i,j) = alpha*b(i,j)
  220                     CONTINUE
                      END IF
                      DO 240 k = j + 1,n
                          IF (a(k,j).NE.zero) THEN
                              DO 230 i = 1,m
                                  b(i,j) = b(i,j) - a(k,j)*b(i,k)
  230                         CONTINUE
                          END IF
  240                 CONTINUE
                      IF (nounit) THEN
!                          temp = one/a(j,j)
                          atmp=0.d0
                          if(a(j,j).gt.tol)atmp=1.d0/a(j,j)
                          temp = one*atmp
                          DO 250 i = 1,m
                              b(i,j) = temp*b(i,j)
  250                     CONTINUE
                      END IF
  260             CONTINUE
              END IF
          ELSE
C
C           Form  B := alpha*B*inv( A**T ).
C
              IF (upper) THEN
                  DO 310 k = n,1,-1
                      IF (nounit) THEN
!                          temp = one/a(k,k)
                          atmp=0.d0
                          if(a(k,k).gt.tol)atmp=1.d0/a(k,k)
                          temp = one*atmp
                          DO 270 i = 1,m
                              b(i,k) = temp*b(i,k)
  270                     CONTINUE
                      END IF
                      DO 290 j = 1,k - 1
                          IF (a(j,k).NE.zero) THEN
                              temp = a(j,k)
                              DO 280 i = 1,m
                                  b(i,j) = b(i,j) - temp*b(i,k)
  280                         CONTINUE
                          END IF
  290                 CONTINUE
                      IF (alpha.NE.one) THEN
                          DO 300 i = 1,m
                              b(i,k) = alpha*b(i,k)
  300                     CONTINUE
                      END IF
  310             CONTINUE
              ELSE
                  DO 360 k = 1,n
                      IF (nounit) THEN
!                          temp = one/a(k,k)
                          atmp=0.d0
                          if(a(k,k).gt.tol)atmp=1.d0/a(k,k)
                          temp = one*atmp
                          DO 320 i = 1,m
                              b(i,k) = temp*b(i,k)
  320                     CONTINUE
                      END IF
                      DO 340 j = k + 1,n
                          IF (a(j,k).NE.zero) THEN
                              temp = a(j,k)
                              DO 330 i = 1,m
                                  b(i,j) = b(i,j) - temp*b(i,k)
  330                         CONTINUE
                          END IF
  340                 CONTINUE
                      IF (alpha.NE.one) THEN
                          DO 350 i = 1,m
                              b(i,k) = alpha*b(i,k)
  350                     CONTINUE
                      END IF
  360             CONTINUE
              END IF
          END IF
      END IF
C
      RETURN
C
C     End of DTRSM .
C
      END

ext: f


text:
*DECK RADB2
      SUBROUTINE RADB2 (IDO, L1, CC, CH, WA1)
C***BEGIN PROLOGUE  RADB2
C***SUBSIDIARY
C***PURPOSE  Calculate the fast Fourier transform of subvectors of
C            length two.
C***LIBRARY   SLATEC (FFTPACK)
C***TYPE      SINGLE PRECISION (RADB2-S)
C***AUTHOR  Swarztrauber, P. N., (NCAR)
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790601  DATE WRITTEN
C   830401  Modified to use SLATEC library source file format.
C   860115  Modified by Ron Boisvert to adhere to Fortran 77 by
C           changing dummy array size declarations (1) to (*).
C   881128  Modified by Dick Valent to meet prologue standards.
C   890831  Modified array declarations.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  RADB2
      DIMENSION CC(IDO,2,*), CH(IDO,L1,2), WA1(*)
C***FIRST EXECUTABLE STATEMENT  RADB2
      DO 101 K=1,L1
         CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
         CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
  101 CONTINUE
      IF (IDO-2) 107,105,102
  102 IDP2 = IDO+2
      IF((IDO-1)/2.LT.L1) GO TO 108
      DO 104 K=1,L1
CDIR$ IVDEP
         DO 103 I=3,IDO,2
            IC = IDP2-I
            CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
            TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
            CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
            TI2 = CC(I,1,K)+CC(IC,2,K)
            CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
            CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
  103    CONTINUE
  104 CONTINUE
      GO TO 111
  108 DO 110 I=3,IDO,2
         IC = IDP2-I
CDIR$ IVDEP
         DO 109 K=1,L1
            CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
            TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
            CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
            TI2 = CC(I,1,K)+CC(IC,2,K)
            CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
            CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
  109    CONTINUE
  110 CONTINUE
  111 IF (MOD(IDO,2) .EQ. 1) RETURN
  105 DO 106 K=1,L1
         CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
         CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
  106 CONTINUE
  107 RETURN
      END

max_stars_repo_name: homiak/jabaws

max_issues_repo_name: homiak/jabaws

max_forks_repo_name: homiak/jabaws


text:
SUBROUTINE classico(uint,vint,wint)

	USE velpre
	USE parametros

	IMPLICIT NONE

	!===================================================================================================================
	real(8), dimension(nx1,ny,nz) :: uint
	real(8), dimension(nx,ny1,nz) :: vint
	real(8), dimension(nx,ny,nz1) :: wint

	real(8), dimension(nx1,ny,nz) :: dudx, dudy, dudz, bma, dma, dudxa, dudya, dudza
	real(8), dimension(nx,ny1,nz) :: dvdx, dvdy, dvdz, amb, dmb, dvdxa, dvdya, dvdza
	real(8), dimension(nx,ny,nz1) :: dwdx, dwdy, dwdz, amd, bmd, dwdxa, dwdya, dwdza
	real(8), dimension(nx,ny,nz)  :: aux

	!contadores
	integer :: i, j, k

	!auxiliares
	real(8) :: aux1, aux2

	!===================================================================================================================
	!RESOLUÇÃO DO PROBLEMA
	!===================================================================================================================

	if (der == 1) then
	! upwind

	do k = 1, nz
	do j = 1, ny
	do i = 1, nx1

	bma(i,j,k) = (v(i,j,k) + v(i-1,j,k) + v(i,j+1,k) + v(i-1,j+1,k)) * 0.25
	dma(i,j,k) = (w(i,j,k) + w(i-1,j,k) + w(i,j,k+1) + w(i-1,j,k+1)) * 0.25

	aux1 = max(u(i,j,k),0.)
	aux2 = min(u(i,j,k),0.)
	dudx(i,j,k) = aux1*(u(i,j,k)-u(i-1,j,k))/dx + aux2*(u(i+1,j,k)-u(i,j,k))/dx

	aux1 = max(bma(i,j,k),0.)
	aux2 = min(bma(i,j,k),0.)
	dudy(i,j,k) = aux1*(u(i,j,k)-u(i,j-1,k))/dy + aux2*(u(i,j+1,k)-u(i,j,k))/dy

	aux1 = max(dma(i,j,k),0.)
	aux2 = min(dma(i,j,k),0.)
	dudz(i,j,k) = aux1*(u(i,j,k)-u(i,j,k-1))/dz + aux2*(u(i,j,k+1)-u(i,j,k))/dz

	uint(i,j,k) = dudx(i,j,k) + dudy(i,j,k) + dudz(i,j,k)

	enddo
	enddo
	enddo

	do k = 1, nz
	do j = 1, ny1
	do i = 1, nx

	amb(i,j,k) = (u(i,j,k) + u(i+1,j,k) + u(i,j-1,k) + u(i+1,j-1,k)) * 0.25
	dmb(i,j,k) = (w(i,j,k) + w(i,j-1,k) + w(i,j,k+1) + w(i,j-1,k+1)) * 0.25

	aux1 = max(amb(i,j,k),0.)
	aux2 = min(amb(i,j,k),0.)
	dvdx(i,j,k) = aux1*(v(i,j,k)-v(i-1,j,k))/dx + aux2*(v(i+1,j,k)-v(i,j,k))/dx

	aux1 = max(v(i,j,k),0.)
	aux2 = min(v(i,j,k),0.)
	dvdy(i,j,k) = aux1*(v(i,j,k)-v(i,j-1,k))/dy + aux2*(v(i,j+1,k)-v(i,j,k))/dy

	aux1 = max(dmb(i,j,k),0.)
	aux2 = min(dmb(i,j,k),0.)
	dvdz(i,j,k) = aux1*(v(i,j,k)-v(i,j,k-1))/dz + aux2*(v(i,j,k+1)-v(i,j,k))/dz

	vint(i,j,k) = dvdx(i,j,k) + dvdy(i,j,k) + dvdz(i,j,k)

	enddo
	enddo
	enddo

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

	amd(i,j,k) = (u(i,j,k) + u(i+1,j,k) + u(i,j,k-1) + u(i+1,j,k-1)) * 0.25
	bmd(i,j,k) = (v(i,j,k) + v(i,j+1,k) + v(i,j,k-1) + v(i,j+1,k-1)) * 0.25

	aux1 = max(amd(i,j,k),0.)
	aux2 = min(amd(i,j,k),0.)
	dwdx(i,j,k) = aux1*(w(i,j,k)-w(i-1,j,k))/dx + aux2*(w(i+1,j,k)-w(i,j,k))/dx

	aux1 = max(bmd(i,j,k),0.)
	aux2 = min(bmd(i,j,k),0.)
	dwdy(i,j,k) = aux1*(w(i,j,k)-w(i,j-1,k))/dy + aux2*(w(i,j+1,k)-w(i,j,k))/dy

	aux1 = max(w(i,j,k),0.)
	aux2 = min(w(i,j,k),0.)
	dwdz(i,j,k) = aux1*(w(i,j,k)-w(i,j,k-1))/dz + aux2*(w(i,j,k+1)-w(i,j,k))/dz

	wint(i,j,k) = dwdx(i,j,k) + dwdy(i,j,k) + dwdz(i,j,k)

	enddo
	enddo
	enddo

	elseif (der == 2) then

	CALL derivax(u,nx1,ny,nz,dudx)
	CALL derivay(u,nx1,ny,nz,dudy)
	CALL derivaz(u,nx1,ny,nz,dudz)

	CALL derivax(v,nx,ny1,nz,dvdx)
	CALL derivay(v,nx,ny1,nz,dvdy)
	CALL derivaz(v,nx,ny1,nz,dvdz)

	CALL derivax(w,nx,ny,nz1,dwdx)
	CALL derivay(w,nx,ny,nz1,dwdy)
	CALL derivaz(w,nx,ny,nz1,dwdz)

	call interpy_fc(v(1:nx,1:ny1,1:nz),nx,ny1,nz,aux) !(nx,ny,nz)
	call interpx_cf(aux,nx,ny,nz,bma) !(nx1,ny,nz)

	call interpz_fc(w(1:nx,1:ny,1:nz1),nx,ny,nz1,aux) !(nx,ny,nz)
	call interpx_cf(aux,nx,ny,nz,dma) !(nx1,ny,nz)

	call interpx_fc(u(1:nx1,1:ny,1:nz),nx1,ny,nz,aux) !(nx,ny,nz)
	call interpy_cf(aux,nx,ny,nz,amb) !(nx,ny1,nz)

	call interpz_fc(w(1:nx,1:ny,1:nz1),nx,ny,nz1,aux) !(nx,ny,nz)
	call interpy_cf(aux,nx,ny,nz,dmb) !(nx,ny1,nz)

	call interpx_fc(u(1:nx1,1:ny,1:nz),nx1,ny,nz,aux) !(nx,ny,nz)
	call interpz_cf(aux,nx,ny,nz,amd) !(nx,ny,nz1)

	call interpy_fc(v(1:nx,1:ny1,1:nz),nx,ny1,nz,aux) !(nx,ny,nz)
	call interpz_cf(aux,nx,ny,nz,bmd) !(nx,ny,nz1)

	do k = 1, nz
	do j = 1, ny
	do i = 1, nx1
	uint(i,j,k) = u(i,j,k)*dudx(i,j,k) + bma(i,j,k)*dudy(i,j,k) + dma(i,j,k)*dudz(i,j,k)
	enddo
	enddo
	enddo

	do k = 1, nz
	do j = 1, ny1
	do i = 1, nx
	vint(i,j,k) = amb(i,j,k)*dvdx(i,j,k) + v(i,j,k)*dvdy(i,j,k) + dmb(i,j,k)*dvdz(i,j,k)
	enddo
	enddo
	enddo

	do k = 1, nz1
	do j = 1, ny
	do i = 1, nx
	wint(i,j,k) = amd(i,j,k)*dwdx(i,j,k) + bmd(i,j,k)*dwdy(i,j,k) + w(i,j,k)*dwdz(i,j,k)
	enddo
	enddo
	enddo

	elseif (der == 3) then
	! upwind 2nd order

	call interpy_fc(v(1:nx,1:ny1,1:nz),nx,ny1,nz,aux) !(nx,ny,nz)
	call interpx_cf(aux,nx,ny,nz,bma) !(nx1,ny,nz)

	call interpz_fc(w(1:nx,1:ny,1:nz1),nx,ny,nz1,aux) !(nx,ny,nz)
	call interpx_cf(aux,nx,ny,nz,dma) !(nx1,ny,nz)

	call interpx_fc(u(1:nx1,1:ny,1:nz),nx1,ny,nz,aux) !(nx,ny,nz)
	call interpy_cf(aux,nx,ny,nz,amb) !(nx,ny1,nz)

	call interpz_fc(w(1:nx,1:ny,1:nz1),nx,ny,nz1,aux) !(nx,ny,nz)
	call interpy_cf(aux,nx,ny,nz,dmb) !(nx,ny1,nz)

	call interpx_fc(u(1:nx1,1:ny,1:nz),nx1,ny,nz,aux) !(nx,ny,nz)
	call interpz_cf(aux,nx,ny,nz,amd) !(nx,ny,nz1)

	call interpy_fc(v(1:nx,1:ny1,1:nz),nx,ny1,nz,aux) !(nx,ny,nz)
	call interpz_cf(aux,nx,ny,nz,bmd) !(nx,ny,nz1)

	CALL derivaxu2n(u,nx1,ny,nz,dudxa)
	CALL derivaxu2p(u,nx1,ny,nz,dudx)
	CALL derivayu2n(u,nx1,ny,nz,dudya)
	CALL derivayu2p(u,nx1,ny,nz,dudy)
	CALL derivazu2n(u,nx1,ny,nz,dudza)
	CALL derivazu2p(u,nx1,ny,nz,dudz)

	CALL derivaxu2n(v,nx,ny1,nz,dvdxa)
	CALL derivaxu2p(v,nx,ny1,nz,dvdx)
	CALL derivayu2n(v,nx,ny1,nz,dvdya)
	CALL derivayu2p(v,nx,ny1,nz,dvdy)
	CALL derivazu2n(v,nx,ny1,nz,dvdza)
	CALL derivazu2p(v,nx,ny1,nz,dvdz)

	CALL derivaxu2n(w,nx,ny,nz1,dwdxa)
	CALL derivaxu2p(w,nx,ny,nz1,dwdx)
	CALL derivayu2n(w,nx,ny,nz1,dwdya)
	CALL derivayu2p(w,nx,ny,nz1,dwdy)
	CALL derivazu2n(w,nx,ny,nz1,dwdza)
	CALL derivazu2p(w,nx,ny,nz1,dwdz)

	do k = 1, nz
	do j = 1, ny
	do i = 1, nx1

	aux1 = max(u(i,j,k),0.)
	aux2 = min(u(i,j,k),0.)
	dudx(i,j,k) = aux1*dudxa(i,j,k) + aux2*dudx(i,j,k)

	aux1 = max(bma(i,j,k),0.)
	aux2 = min(bma(i,j,k),0.)
	dudy(i,j,k) = aux1*dudya(i,j,k) + aux2*dudy(i,j,k)

	aux1 = max(dma(i,j,k),0.)
	aux2 = min(dma(i,j,k),0.)
	dudz(i,j,k) = aux1*dudza(i,j,k) + aux2*dudz(i,j,k)

	uint(i,j,k) = dudx(i,j,k) + dudy(i,j,k) + dudz(i,j,k)

	enddo
	enddo
	enddo

	do k = 1, nz
	do j = 1, ny1
	do i = 1, nx

	aux1 = max(amb(i,j,k),0.)
	aux2 = min(amb(i,j,k),0.)
	dvdx(i,j,k) = aux1*dvdxa(i,j,k) + aux2*dvdx(i,j,k)

	aux1 = max(v(i,j,k),0.)
	aux2 = min(v(i,j,k),0.)
	dvdy(i,j,k) = aux1*dvdya(i,j,k) + aux2*dvdy(i,j,k)

	aux1 = max(dmb(i,j,k),0.)
	aux2 = min(dmb(i,j,k),0.)
	dvdz(i,j,k) = aux1*dvdza(i,j,k) + aux2*dvdz(i,j,k)

	vint(i,j,k) = dvdx(i,j,k) + dvdy(i,j,k) + dvdz(i,j,k)

	enddo
	enddo
	enddo

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

	aux1 = max(amd(i,j,k),0.)
	aux2 = min(amd(i,j,k),0.)
	dwdx(i,j,k) = aux1*dwdxa(i,j,k) + aux2*dwdx(i,j,k)

	aux1 = max(bmd(i,j,k),0.)
	aux2 = min(bmd(i,j,k),0.)
	dwdy(i,j,k) = aux1*dwdya(i,j,k) + aux2*dwdy(i,j,k)

	aux1 = max(w(i,j,k),0.)
	aux2 = min(w(i,j,k),0.)
	dwdz(i,j,k) = aux1*dwdza(i,j,k) + aux2*dwdz(i,j,k)

	wint(i,j,k) = dwdx(i,j,k) + dwdy(i,j,k) + dwdz(i,j,k)

	enddo
	enddo
	enddo

	endif

!==================================================================================================================
END SUBROUTINE classico

SUBROUTINE rotacional(uint,vint,wint)

	USE velpre
	USE parametros

	IMPLICIT NONE

	!===================================================================================================================
	real(8), dimension(nx1,ny,nz) :: uint
	real(8), dimension(nx,ny1,nz) :: vint
	real(8), dimension(nx,ny,nz1) :: wint

	real(8), dimension(0:nx1+1,0:ny1+1,0:nz1+1) :: ap, an, bma, dma
	real(8), dimension(0:nx1+1,0:ny1+1,0:nz1+1) :: bp, bn, amb, dmb
	real(8), dimension(0:nx1+1,0:ny1+1,0:nz1+1) :: dp, dn, amd, bmd

	!
	real(8), dimension(nx1,ny1,nz1) :: dudx, dvdx, dwdx
	real(8), dimension(nx1,ny1,nz1) :: dudy, dvdy, dwdy
	real(8), dimension(nx1,ny1,nz1) :: dudz, dvdz, dwdz

	!
	real(8) :: aa, bb, dd

	!contadores
	integer :: i, j, k, ai, bi, di

	!plotagem
	real(8) :: acont, bcont, dcont
	integer :: loca(3), locb(3), locd(3)

	!auxiliares
	real(8) :: aux1, aux2

	!===================================================================================================================
	!RESOLUÇÃO DO PROBLEMA
	!===================================================================================================================

	if (der == 1) then
	! upwind

	do k = 0, nz+1
	do j = 0, ny+1
	do i = 1, nx+1

	bma(i,j,k) = (v(i,j,k) + v(i-1,j,k) + v(i,j+1,k) + v(i-1,j+1,k)) * 0.25
	dma(i,j,k) = (w(i,j,k) + w(i-1,j,k) + w(i,j,k+1) + w(i-1,j,k+1)) * 0.25

	ap(i,j,k) = max(u(i,j,k),0.)
	an(i,j,k) = min(u(i,j,k),0.)

	bp(i,j,k) = max(bma(i,j,k),0.)
	bn(i,j,k) = min(bma(i,j,k),0.)

	dp(i,j,k) = max(dma(i,j,k),0.)
	dn(i,j,k) = min(dma(i,j,k),0.)

	enddo
	enddo
	enddo

	ap(0,0:ny+1,0:nz+1)     = max(u(0,0:ny+1,0:nz+1),0.)
	an(nx1+1,0:ny+1,0:nz+1) = min(u(nx1+1,0:ny+1,0:nz+1),0.)

	do k = 1, nz
	do j = 1, ny
	do i = 1, nx1

	dudx(i,j,k) = (ap(i,j,k)*u(i,j,k)-ap(i-1,j,k)*u(i-1,j,k))/dx + (an(i+1,j,k)*u(i+1,j,k)-an(i,j,k)*u(i,j,k))/dx
	dudy(i,j,k) = (bp(i,j,k)*u(i,j,k)-bp(i,j-1,k)*u(i,j-1,k))/dy + (bn(i,j+1,k)*u(i,j+1,k)-bn(i,j,k)*u(i,j,k))/dy
	dudz(i,j,k) = (dp(i,j,k)*u(i,j,k)-dp(i,j,k-1)*u(i,j,k-1))/dz + (dn(i,j,k+1)*u(i,j,k+1)-dn(i,j,k)*u(i,j,k))/dz

	uint(i,j,k) = dudx(i,j,k) + dudy(i,j,k) + dudz(i,j,k)

	enddo
	enddo
	enddo

	do k = 0, nz+1
	do j = 1, ny+1
	do i = 0, nx+1

	amb(i,j,k) = (u(i,j,k) + u(i+1,j,k) + u(i,j-1,k) + u(i+1,j-1,k)) * 0.25
	dmb(i,j,k) = (w(i,j,k) + w(i,j-1,k) + w(i,j,k+1) + w(i,j-1,k+1)) * 0.25

	ap(i,j,k) = max(amb(i,j,k),0.)
	an(i,j,k) = min(amb(i,j,k),0.)

	bp(i,j,k) = max(v(i,j,k),0.)
	bn(i,j,k) = min(v(i,j,k),0.)

	dp(i,j,k) = max(dmb(i,j,k),0.)
	dn(i,j,k) = min(dmb(i,j,k),0.)

	enddo
	enddo
	enddo

	bp(0:nx+1,0,0:nz+1)     = max(v(0:nx+1,0,0:nz+1),0.)
	bn(0:nx+1,ny1+1,0:nz+1) = min(v(0:nx+1,ny1+1,0:nz+1),0.)

	do k = 1, nz
	do j = 1, ny1
	do i = 1, nx

	dvdx(i,j,k) = (ap(i,j,k)*v(i,j,k)-ap(i-1,j,k)*v(i-1,j,k))/dx + (an(i+1,j,k)*v(i+1,j,k)-an(i,j,k)*v(i,j,k))/dx
	dvdy(i,j,k) = (bp(i,j,k)*v(i,j,k)-bp(i,j-1,k)*v(i,j-1,k))/dy + (bn(i,j+1,k)*v(i,j+1,k)-bn(i,j,k)*v(i,j,k))/dy
	dvdz(i,j,k) = (dp(i,j,k)*v(i,j,k)-dp(i,j,k-1)*v(i,j,k-1))/dz + (dn(i,j,k+1)*v(i,j,k+1)-dn(i,j,k)*v(i,j,k))/dz

	vint(i,j,k) = dvdx(i,j,k) + dvdy(i,j,k) + dvdz(i,j,k)

	enddo
	enddo
	enddo

	do k = 1, nz+1
	do j = 0, ny+1
	do i = 0, nx+1

	amd(i,j,k) = (u(i,j,k) + u(i+1,j,k) + u(i,j,k-1) + u(i+1,j,k-1)) * 0.25
	bmd(i,j,k) = (v(i,j,k) + v(i,j+1,k) + v(i,j,k-1) + v(i,j+1,k-1)) * 0.25

	ap(i,j,k) = max(amd(i,j,k),0.)
	an(i,j,k) = min(amd(i,j,k),0.)

	bp(i,j,k) = max(bmd(i,j,k),0.)
	bn(i,j,k) = min(bmd(i,j,k),0.)

	dp(i,j,k) = max(w(i,j,k),0.)
	dn(i,j,k) = min(w(i,j,k),0.)

	enddo
	enddo
	enddo

	dp(0:nx+1,0:ny+1,0)     = max(w(0:nx+1,0:ny+1,0),0.)
	dp(0:nx+1,0:ny+1,nz1+1) = min(w(0:nx+1,0:ny+1,nz1+1),0.)
	do k = 1, nz1
	do j = 1, ny
	do i = 1, nx

	dwdx(i,j,k) = (ap(i,j,k)*w(i,j,k)-ap(i-1,j,k)*w(i-1,j,k))/dx + (an(i+1,j,k)*w(i+1,j,k)-an(i,j,k)*w(i,j,k))/dx
	dwdy(i,j,k) = (bp(i,j,k)*w(i,j,k)-bp(i,j-1,k)*w(i,j-1,k))/dy + (bn(i,j+1,k)*w(i,j+1,k)-bn(i,j,k)*w(i,j,k))/dy
	dwdz(i,j,k) = (dp(i,j,k)*w(i,j,k)-dp(i,j,k-1)*w(i,j,k-1))/dz + (dn(i,j,k+1)*w(i,j,k+1)-dn(i,j,k)*w(i,j,k))/dz

	wint(i,j,k) = dwdx(i,j,k) + dwdy(i,j,k) + dwdz(i,j,k)

	enddo
	enddo
	enddo

	else

	write(*,*) "não possui este esquema para o rotacional"
	STOP

	endif

!==================================================================================================================
END SUBROUTINE rotacional

SUBROUTINE antissim(uint,vint,wint)

	USE velpre
	USE parametros

	IMPLICIT NONE

	!===================================================================================================================
	real(8), dimension(nx1,ny,nz) :: uint
	real(8), dimension(nx,ny1,nz) :: vint
	real(8), dimension(nx,ny,nz1) :: wint

	real(8), dimension(0:nx1+1,0:ny1+1,0:nz1+1) :: ap, an
	real(8), dimension(0:nx1+1,0:ny1+1,0:nz1+1) :: bp, bn
	real(8), dimension(0:nx1+1,0:ny1+1,0:nz1+1) :: dp, dn

	!
	real(8), dimension(nx1,ny1,nz1) :: dudx, dvdx, dwdx
	real(8), dimension(nx1,ny1,nz1) :: dudy, dvdy, dwdy
	real(8), dimension(nx1,ny1,nz1) :: dudz, dvdz, dwdz
	real(8), dimension(nx1,ny1,nz1) :: bma, dma,amb, dmb, amd, bmd

	!
	real(8) :: aa, bb, dd

	!contadores
	integer :: i, j, k, ai, bi, di

	!plotagem
	real(8) :: acont, bcont, dcont
	integer :: loca(3), locb(3), locd(3)

	!auxiliares
	real(8) :: aux1, aux2

	!===================================================================================================================
	!RESOLUÇÃO DO PROBLEMA
	!===================================================================================================================

	if (der == 1) then
	! upwind

	do k = 1, nz+1
	do j = 1, ny+1
	do i = 1, nx+1
	bma(i,j,k) = (v(i,j,k) + v(i-1,j,k) + v(i,j+1,k) + v(i-1,j+1,k)) * 0.25
	dma(i,j,k) = (w(i,j,k) + w(i-1,j,k) + w(i,j,k+1) + w(i-1,j,k+1)) * 0.25

	ap(i,j,k) = max(u(i,j,k),0.)
	an(i,j,k) = min(u(i,j,k),0.)

	bp(i,j,k) = max(bma(i,j,k),0.)
	bn(i,j,k) = min(bma(i,j,k),0.)

	dp(i,j,k) = max(dma(i,j,k),0.)
	dn(i,j,k) = min(dma(i,j,k),0.)
	enddo
	enddo
	enddo

	ap(0,:,:) = ap(1,:,:)
	bp(:,0,:) = bp(:,1,:)
	dp(:,:,0) = dp(:,:,1)
        an(nx1+1,:,:)=an(nx1,:,:)

	do k = 1, nz
	do j = 1, ny
	do i = 1, nx1

	dudx(i,j,k) = ((ap(i,j,k)*u(i,j,k)-ap(i-1,j,k)*u(i-1,j,k))/dx + (an(i+1,j,k)*u(i+1,j,k)-an(i,j,k)*u(i,j,k))/dx &
	+ ap(i,j,k)*(u(i,j,k)-u(i-1,j,k))/dx + an(i,j,k)*(u(i+1,j,k)-u(i,j,k))/dx) * 0.5

	dudy(i,j,k) = ((bp(i,j,k)*u(i,j,k)-bp(i,j-1,k)*u(i,j-1,k))/dy + (bn(i,j+1,k)*u(i,j+1,k)-bn(i,j,k)*u(i,j,k))/dy &
	+ bp(i,j,k)*(u(i,j,k)-u(i,j-1,k))/dy + bn(i,j,k)*(u(i,j+1,k)-u(i,j,k))/dy) * 0.5

	dudz(i,j,k) = ((dp(i,j,k)*u(i,j,k)-dp(i,j,k-1)*u(i,j,k-1))/dz+(dn(i,j,k+1)*u(i,j,k+1)-dn(i,j,k)*u(i,j,k))/dz &
	+ dp(i,j,k)*(u(i,j,k)-u(i,j,k-1))/dz + dn(i,j,k)*(u(i,j,k+1)-u(i,j,k))/dz) * 0.5

	uint(i,j,k) = dudx(i,j,k) + dudy(i,j,k) + dudz(i,j,k)

	enddo
	enddo
	enddo

	do k = 1, nz+1
	do j = 1, ny+1
	do i = 1, nx+1

	amb(i,j,k) = (u(i,j,k) + u(i+1,j,k) + u(i,j-1,k) + u(i+1,j-1,k)) * 0.25
	dmb(i,j,k) = (w(i,j,k) + w(i,j-1,k) + w(i,j,k+1) + w(i,j-1,k+1)) * 0.25

	ap(i,j,k) = max(amb(i,j,k),0.)
	an(i,j,k) = min(amb(i,j,k),0.)

	bp(i,j,k) = max(v(i,j,k),0.)
	bn(i,j,k) = min(v(i,j,k),0.)

	dp(i,j,k) = max(dmb(i,j,k),0.)
	dn(i,j,k) = min(dmb(i,j,k),0.)

	enddo
	enddo
	enddo

	ap(0,:,:) = ap(1,:,:)
	bp(:,0,:) = bp(:,1,:)
	dp(:,:,0) = dp(:,:,1)
        bn(:,ny1+1,:)=bn(:,ny1,:)

	do k = 1, nz
	do j = 1, ny1
	do i = 1, nx

	dvdx(i,j,k) = ((ap(i,j,k)*v(i,j,k)-ap(i-1,j,k)*v(i-1,j,k))/dx + (an(i+1,j,k)*v(i+1,j,k)-an(i,j,k)*v(i,j,k))/dx &
	+ ap(i,j,k)*(v(i,j,k)-v(i-1,j,k))/dx + an(i,j,k)*(v(i+1,j,k)-v(i,j,k))/dx) * 0.5

	dvdy(i,j,k) = ((bp(i,j,k)*v(i,j,k)-bp(i,j-1,k)*v(i,j-1,k))/dy + (bn(i,j+1,k)*v(i,j+1,k)-bn(i,j,k)*v(i,j,k))/dy &
	+ bp(i,j,k)*(v(i,j,k)-v(i,j-1,k))/dy + bn(i,j,k)*(v(i,j+1,k)-v(i,j,k))/dy) * 0.5

	dvdz(i,j,k) = ((dp(i,j,k)*v(i,j,k)-dp(i,j,k-1)*v(i,j,k-1))/dz+(dn(i,j,k+1)*v(i,j,k+1)-dn(i,j,k)*v(i,j,k))/dz &
	+ dp(i,j,k)*(v(i,j,k)-v(i,j,k-1))/dz + dn(i,j,k)*(v(i,j,k+1)-v(i,j,k))/dz) * 0.5

	vint(i,j,k) = dvdx(i,j,k) + dvdy(i,j,k) + dvdz(i,j,k)

	enddo
	enddo
	enddo

	do k = 1, nz+1
	do j = 1, ny+1
	do i = 1, nx+1

	amd(i,j,k) = (u(i,j,k) + u(i+1,j,k) + u(i,j,k-1) + u(i+1,j,k-1)) * 0.25
	bmd(i,j,k) = (v(i,j,k) + v(i,j+1,k) + v(i,j,k-1) + v(i,j+1,k-1)) * 0.25

	ap(i,j,k) = max(amd(i,j,k),0.)
	an(i,j,k) = min(amd(i,j,k),0.)

	bp(i,j,k) = max(bmd(i,j,k),0.)
	bn(i,j,k) = min(bmd(i,j,k),0.)

	dp(i,j,k) = max(w(i,j,k),0.)
	dn(i,j,k) = min(w(i,j,k),0.)

	enddo
	enddo
	enddo

	ap(0,:,:) = ap(1,:,:)
	bp(:,0,:) = bp(:,1,:)
	dp(:,:,0) = dp(:,:,1)
        dn(:,:,nz1+1)=dn(:,:,nz1)

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

	dwdx(i,j,k) = ((ap(i,j,k)*w(i,j,k)-ap(i-1,j,k)*w(i-1,j,k))/dx + (an(i+1,j,k)*w(i+1,j,k)-an(i,j,k)*w(i,j,k))/dx &
	+ ap(i,j,k)*(w(i,j,k)-w(i-1,j,k))/dx + an(i,j,k)*(w(i+1,j,k)-w(i,j,k))/dx) * 0.5

	dwdy(i,j,k) = ((bp(i,j,k)*w(i,j,k)-bp(i,j-1,k)*w(i,j-1,k))/dy + (bn(i,j+1,k)*w(i,j+1,k)-bn(i,j,k)*w(i,j,k))/dy &
	+ bp(i,j,k)*(w(i,j,k)-w(i,j-1,k))/dy + bn(i,j,k)*(w(i,j+1,k)-w(i,j,k))/dy) * 0.5

	dwdz(i,j,k) = ((dp(i,j,k)*w(i,j,k)-dp(i,j,k-1)*w(i,j,k-1))/dz+(dn(i,j,k+1)*w(i,j,k+1)-dn(i,j,k)*w(i,j,k))/dz &
	+ dp(i,j,k)*(w(i,j,k)-w(i,j,k-1))/dz + dn(i,j,k)*(w(i,j,k+1)-w(i,j,k))/dz) * 0.5

	wint(i,j,k) = dwdx(i,j,k) + dwdy(i,j,k) + dwdz(i,j,k)

	enddo
	enddo
	enddo

	else

	write(*,*) "não possui esquema para o antissimétrico"
	STOP

	endif

!==================================================================================================================
END SUBROUTINE antissim

max_stars_repo_name: pemiguell/SuLi

max_issues_repo_name: pemiguell/SuLi

max_forks_repo_name: pemiguell/SuLi


text:
	PROGRAM GDCFIL
C************************************************************************
C* GDCFIL								*
C*									*
C* This program creates a new grid file.				*
C*									*
C* Log: 								*
C* I. Graffman/RDS	 6/88						*
C* M. desJardins/GSFC	 8/88	Lots of changes				*
C* G. Huffman/GSC	 4/89	Grid navigation table input		*
C* M. desJardins/GSFC	 5/89	Error messages; init ihd, anlblk	*
C* K. Brill/GSC          4/90   Added ANLYSS for ANLBLK input		*
C* K. Brill/NMC          7/90   Changes for table file listing		*
C* S. Schotz/GSC	10/90	Eliminated invalid error message	*
C* M. desJardins/NMC	 4/91	Eliminated xspace, yspace		*
C* J. Whistler/SSAI	 5/91	Added error message for large grids	*
C* K. Brill/NMC		02/92	Use LLNANL, LLNNAV			*
C* K. Brill/NMC		02/92	Use CPYFIL analysis block only when no	*
C*			   	analysis block is specified.		*
C* L. Williams/EAi      03/94   Clean up declarations of user input	*
C*				variables				*
C* L. Williams/EAI	 7/94	Removed call to GDCUPD			*
C* K. Tyle/GSC		 8/96	Added FL_MFIL to search for file type	*
C* S. Maxwell/GSC	 7/97	Increased input character length 	*
C* T. Lee/GSC		 7/99	Checked CPYFIL				*
C* D.W.Plummer/NCEP	 3/00	Added subarea option to CPYFIL		*
C* S. Jacobs/NCEP	 3/01	Changed FL_MFIL to FL_MFLS		*
C* K. Brill/HPC		11/02	Check KX * KY against LLMXTG		*
C* M. Li/SAIC		04/04	set ihd = LLGDHD			*
C* R. Tian/SAIC          3/05   Changes for time/file mngmnt            *
C* T. Piper/SAIC        01/08   Added GD_INIT; removed from IN_BDTA     *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER	gdfile*(LLMXLN), proj*(LLMXLN), cpyfil*(LLMXLN),
     +			gdarea*(LLMXLN), anlyss*(LLMXLN), kxky*(LLMXLN),
     +		 	maxgrd*(LLMXLN)
C*
	LOGICAL		respnd, done, proces
C-----------------------------------------------------------------------
C*	Initialize TAE.
C
	CALL IP_INIT  ( respnd, iperr )
	IF  ( iperr .eq. 0 )  THEN
	    CALL IP_IDNT  ( 'GDCFIL', ier )
	    done = .false.
	  ELSE
	    done = .true.
	END IF
	IF  ( .not. done )  THEN
C
C*      Initialize grid library common area grdcmn.cmn
C
            CALL GD_INIT  ( iret )
C
C*      Initialize GEMPLT in order to set grid navigation later.
C
	    CALL GG_INIT  ( 1, iret )
	    IF  ( iret .ne. 0 )  THEN
		iperr = -3
		done  = .true.
	    END IF
	END IF
C
C*	Create new grid files.
C
	DO WHILE  ( .not. done )
C
C*	    Get user input.
C
	    CALL GDCINP  ( gdfile, proj, gdarea, kxky, maxgrd, cpyfil,
     +			   anlyss, iperr )
	    IF  ( iperr .ne. 0 )  THEN
	        done   = .true.
		proces = .false.
	      ELSE
	        proces = .true.
	    END IF
C
C*	    Check that output file name is not blank.
C
	    IF  ( proces )  THEN
		IF  ( gdfile .eq. ' ' )  THEN
		    iret = -7
		    CALL ER_WMSG  ( 'GDCFIL', iret, gdfile, ier )
		    proces = .false.
		END IF
	    END IF
C
C*	    Give user option to exit.
C
	    IF  ( proces )  THEN
		CALL GDCDSP  ( gdfile, proj, gdarea, kxky, maxgrd,
     +			       cpyfil, anlyss, iret )
		IF  ( iret .ne. 0 )  proces = .false.
	    END IF
C
C*	    Create the grid file.
C
	    IF  ( proces )  THEN
		CALL GDGCFL  ( gdfile, proj, gdarea, kxky, maxgrd,
     +				cpyfil, anlyss, iret )
	        IF ( iret .ne. 0 ) THEN
		    CALL ER_WMSG  ( 'GDCFIL', iret, gdfile, ier )
		END IF
	    END IF
C
C*	    Call dynamic tutor.
C
	    CALL IP_DYNM  ( done, ier )
	END DO
C
C*	Final error messages.
C
	IF  ( iperr .ne. 0 )  CALL ER_WMSG ( 'GDCFIL', iperr, ' ', ier )
	CALL IP_EXIT  ( iret )
C*
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
module linear_solve

   implicit none

   public :: lu_decomposition
   public :: lu_back_substitution
   public :: lu_inverse
   public :: imaxloc

   interface lu_decomposition
      module procedure lu_decomposition_real
      module procedure lu_decomposition_complex
   end interface

   interface lu_back_substitution
      module procedure lu_back_substitution_real
      module procedure lu_back_substitution_real_complex
      module procedure lu_back_substitution_complex
      module procedure lu_back_substitution_matrix_real
      module procedure lu_back_substitution_matrix_complex
   end interface

   interface lu_inverse
      module procedure lu_inverse_real
      module procedure lu_inverse_complex
   end interface

contains

   subroutine lu_decomposition_real(lu, idx, d)

      implicit none

      real, dimension(:, :), intent(in out) :: lu
      integer, dimension(:), intent(out) :: idx
      real, intent(out) :: d

      real, parameter :: zero = 1.0e-20
      real, dimension(size(lu, 1)) :: vv
      real, dimension(size(lu, 2)) :: dum

      integer :: j, n, imax

      n = size(lu, 1)

      d = 1.0
      vv = maxval(abs(lu), dim=2)
      if (any(vv == 0.0)) &
         write (*, *) 'singular matrix in lu_decomposition'
      vv = 1.0 / vv
      do j = 1, n
         imax = (j - 1) + imaxloc(vv(j:n) * abs(lu(j:n, j)))
         if (j /= imax) then
            dum = lu(imax, :)
            lu(imax, :) = lu(j, :)
            lu(j, :) = dum
            d = -d
            vv(imax) = vv(j)
         end if
         idx(j) = imax
         if (lu(j, j) == 0.0) lu(j, j) = zero
         lu(j + 1:n, j) = lu(j + 1:n, j) / lu(j, j)
         lu(j + 1:n, j + 1:n) = lu(j + 1:n, j + 1:n) - spread(lu(j + 1:n, j), 2, n - j) &
                                * spread(lu(j, j + 1:n), 1, n - j)
      end do

   end subroutine lu_decomposition_real

   subroutine lu_decomposition_complex(lu, idx, d)

      implicit none

      complex, dimension(:, :), intent(in out) :: lu
      integer, dimension(:), intent(out) :: idx
      real, intent(out) :: d

      real, parameter :: zero = 1.0e-20
      real, dimension(size(lu, 1)) :: vv
      complex, dimension(size(lu, 2)) :: dum

      integer :: j, n, imax

      n = size(lu, 1)

      d = 1.0
      vv = maxval(cabs(lu), dim=2)
      if (any(vv == 0.0)) &
         write (*, *) 'singular matrix in lu_decomposition'
      vv = 1.0 / vv
      do j = 1, n
         imax = (j - 1) + imaxloc(vv(j:n) * cabs(lu(j:n, j)))
         if (j /= imax) then
            dum = lu(imax, :)
            lu(imax, :) = lu(j, :)
            lu(j, :) = dum
            d = -d
            vv(imax) = vv(j)
         end if
         idx(j) = imax
         if (lu(j, j) == 0.0) lu(j, j) = zero
         lu(j + 1:n, j) = lu(j + 1:n, j) / lu(j, j)
         lu(j + 1:n, j + 1:n) = lu(j + 1:n, j + 1:n) - spread(lu(j + 1:n, j), 2, n - j) &
                                * spread(lu(j, j + 1:n), 1, n - j)
      end do

   end subroutine lu_decomposition_complex

   subroutine lu_back_substitution_real(lu, idx, b)

      implicit none

      real, dimension(:, :), intent(in) :: lu
      integer, dimension(:), intent(in) :: idx
      real, dimension(:), intent(in out) :: b

      integer :: i, n, ii, ll
      real :: summ

      n = size(lu, 1)
      ii = 0
      do i = 1, n
         ll = idx(i)
         summ = b(ll)
         b(ll) = b(i)
         if (ii /= 0) then
            summ = summ - dot_product(lu(i, ii:i - 1), b(ii:i - 1))
         else if (summ /= 0.0) then
            ii = i
         end if
         b(i) = summ
      end do
      do i = n, 1, -1
         b(i) = (b(i) - dot_product(lu(i, i + 1:n), b(i + 1:n))) / lu(i, i)
      end do

   end subroutine lu_back_substitution_real

   subroutine lu_back_substitution_real_complex(lu, idx, b)

      implicit none

      real, dimension(:, :), intent(in) :: lu
      integer, dimension(:), intent(in) :: idx
      complex, dimension(:), intent(in out) :: b

      integer :: i, n, ii, ll
      complex :: summ

    !! The dot products we use below automatically take the complex conjugate of its
    !! first argument, which we do not want here. Hence the use of conjg to undo this

      n = size(lu, 1)
      ii = 0
      do i = 1, n
         ll = idx(i)
         summ = b(ll)
         b(ll) = b(i)
         if (ii /= 0) then
            summ = summ - dot_product(lu(i, ii:i - 1), b(ii:i - 1))
         else if (summ /= 0.0) then
            ii = i
         end if
         b(i) = summ
      end do
      do i = n, 1, -1
         b(i) = (b(i) - dot_product(lu(i, i + 1:n), b(i + 1:n))) / lu(i, i)
      end do

   end subroutine lu_back_substitution_real_complex

   subroutine lu_back_substitution_complex(lu, idx, b)

      implicit none

      complex, dimension(:, :), intent(in) :: lu
      integer, dimension(:), intent(in) :: idx
      complex, dimension(:), intent(in out) :: b

      integer :: i, n, ii, ll
      complex :: summ

    !! The dot products we use below automatically take the complex conjugate of its
    !! first argument, which we do not want here. Hence the use of conjg to undo this

      n = size(lu, 1)
      ii = 0
      do i = 1, n
         ll = idx(i)
         summ = b(ll)
         b(ll) = b(i)
         if (ii /= 0) then
            summ = summ - dot_product(conjg(lu(i, ii:i - 1)), b(ii:i - 1))
         else if (summ /= 0.0) then
            ii = i
         end if
         b(i) = summ
      end do
      do i = n, 1, -1
         b(i) = (b(i) - dot_product(conjg(lu(i, i + 1:n)), b(i + 1:n))) / lu(i, i)
      end do

   end subroutine lu_back_substitution_complex

   subroutine lu_back_substitution_matrix_real(lu, idx, b)

      implicit none

      real, dimension(:, :), intent(in) :: lu
      integer, dimension(:), intent(in) :: idx
      real, dimension(:, :), intent(in out) :: b

      integer :: i, j, n, ii, ll
      real :: summ

      n = size(lu, 1)
      ii = 0
      do j = 1, n
         do i = 1, n
            ll = idx(i)
            summ = b(ll, j)
            b(ll, j) = b(i, j)
            if (ii /= 0) then
               summ = summ - dot_product(lu(i, ii:i - 1), b(ii:i - 1, j))
            else if (summ /= 0.0) then
               ii = i
            end if
            b(i, j) = summ
         end do
         do i = n, 1, -1
            b(i, j) = (b(i, j) - dot_product(lu(i, i + 1:n), b(i + 1:n, j))) / lu(i, i)
         end do
      end do

   end subroutine lu_back_substitution_matrix_real

   subroutine lu_back_substitution_matrix_complex(lu, idx, b)

      implicit none

      complex, dimension(:, :), intent(in) :: lu
      integer, dimension(:), intent(in) :: idx
      complex, dimension(:, :), intent(in out) :: b

      integer :: i, j, n, ii, ll
      complex :: summ

    !! The dot products we use below automatically take the complex conjugate of its
    !! first argument, which we do not want here. Hence the use of conjg to undo this

      n = size(lu, 1)
      ii = 0
      do j = 1, n
         do i = 1, n
            ll = idx(i)
            summ = b(ll, j)
            b(ll, j) = b(i, j)
            if (ii /= 0) then
               summ = summ - dot_product(conjg(lu(i, ii:i - 1)), b(ii:i - 1, j))
            else if (summ /= 0.0) then
               ii = i
            end if
            b(i, j) = summ
         end do
         do i = n, 1, -1
            b(i, j) = (b(i, j) - dot_product(conjg(lu(i, i + 1:n)), b(i + 1:n, j))) / lu(i, i)
         end do
      end do

   end subroutine lu_back_substitution_matrix_complex

  !!
  !! One shouldn't need to compute the inverse of a matrix if solving the linear equation A.x = y;
  !! a simple back substitution should suffice.
  !!
  !! Only compute the inverse of A when absolutely necessary!
  !!
   subroutine lu_inverse_real(lu, idx, inverse)

      implicit none

      real, dimension(:, :), intent(in) :: lu
      integer, dimension(:), intent(in) :: idx
      real, dimension(:, :), intent(out) :: inverse

      integer :: i, n

      n = size(lu, 1)

      inverse = 0.0
      do i = 1, n
         inverse(i, i) = 1.0
      end do

      do i = 1, n
         call lu_back_substitution(lu, idx, inverse(:, i))
      end do

   end subroutine lu_inverse_real

   subroutine lu_inverse_complex(lu, idx, inverse)

      implicit none

      complex, dimension(:, :), intent(in) :: lu
      integer, dimension(:), intent(in) :: idx
      complex, dimension(:, :), intent(out) :: inverse

      integer :: i, n

      n = size(lu, 1)

      inverse = 0.0
      do i = 1, n
         inverse(i, i) = 1.0
      end do

      do i = 1, n
         call lu_back_substitution(lu, idx, inverse(:, i))
      end do

   end subroutine lu_inverse_complex

   function imaxloc(array)
      real, dimension(:), intent(in) :: array
      integer :: imaxloc
      integer, dimension(1) :: imax
      imax = maxloc(array)
      imaxloc = imax(1)
   end function imaxloc

end module linear_solve

max_stars_repo_name: AntonioG-Jerez/stella

max_issues_repo_name: AntonioG-Jerez/stella

max_forks_repo_name: AntonioG-Jerez/stella


text:
!***********************************************************************
!***********************************************************************
!***********************************************************************
!**                                                                  ***
!**               *******  *******   **    **  *******               ***
!**               **       **    **  **    **  **                    ***
!**               **       **    **  **    **  **                    ***
!**               *****    *******   ** ** **  *****                 ***
!**               **       **  **    ** ** **  **                    ***
!**               **       **   **   ** ** **  **                    ***
!**               *******  **    **   **  **   **                    ***
!**                                                                  ***
!**   Program to create a  .rwf  file for RSCF92 by merging one or   ***
!**   more existing  .rwf  files or generating Thomas-Fermi or hy-   ***
!**   drogenic functions.                                            ***
!**                                                                  ***
!**                              GRASP92                             ***
!**           F. A. Parpia, C. F. Fischer, and I. P. Grant           ***
!**                                                                  ***
!**                                                                  ***
!***********************************************************************
!***********************************************************************
!***********************************************************************
!                                                                      *
      PROGRAM ERWF
!                                                                      *
!   Entry routine for RCI92. Controls the entire computation.          *
!                                                                      *
!   Call(s) to: [LIB92]: SETMC, SETCON.                                *
!               [RCI92]: CHKPLT, MATRIX, SETCSL, SETDBG, SETMIX,       *
!                        SETRES, SETSUM, STRSUM.                       *
!               [NJGRAF]: FACTT.                                       *
!                                                                      *
!   Written by Farid A. Parpia            Last revision: 15 Dec 1992   *
!                                                                      *
!***********************************************************************
!
!...Translated by Pacific-Sierra Research 77to90  4.3E  11:06:21   1/ 2/07
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE DEFAULT_C
      USE CONS_C
      USE IOUNIT_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE getyn_I
      USE setdbg_I
      USE setmc_I
      USE setcon_I
      USE setsum_I
      USE setcsh_I
      USE screenpar_I
      USE getinf_I
      USE strsum_I
      USE genrwf_I
      USE orthsc_I
      USE wrtrwf_I
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: NCORE
      LOGICAL :: YES
!-----------------------------------------------

!   Startup message
!
      WRITE (ISTDE, *) 'RWFNESTIMATE'
      WRITE (ISTDE, *) 'This program estimates radial wave functions'
      WRITE (ISTDE, *) 'for orbitals'
      WRITE (ISTDE, *) 'Input files: isodata, rcsf.inp, optional rwfn file'
      WRITE (ISTDE, *) 'Output file: rwfn.inp'

!
      WRITE (ISTDE, *) 'Default settings ?'
      YES = GETYN()
      IF (YES) THEN
         NDEF = 0
      ELSE
         NDEF = 1
      ENDIF

!
!   Determine if there is to be any debug printout; this will be
!   made on the  .dbg  file
!
      CALL SETDBG
!
!   Perform machine- and installation-dependent setup
!
      CALL SETMC
!
!   Set up the physical constants
!
      CALL SETCON
!
!   Open the  .sum  file
!
      IF (NDEF /= 0) CALL SETSUM
!
!   Open, check, load data from, and close the  .csl  file
!
      CALL SETCSH (21, 'rcsf.inp', NCORE)
!
!   Hydrogenic screen parameters for all orbitals
!
      CALL SCREENPAR (NCORE)
!
!   Determine other relevant information
!
      CALL GETINF
!
!   Write the first part of the  .sum  file
!
      IF (NDEF /= 0) CALL STRSUM
!
!   Generate the subshell radial wavefunctions
!
      CALL GENRWF
!
!   Orthogonalize the radial orbitals
!
      CALL ORTHSC
!
!   Write the subshell radial wavefunctions out
!
      CALL WRTRWF
!
!   Print completion message
!
      WRITE (ISTDE, *) 'RWFNESTIMATE: Execution complete.'
!
      STOP
      END PROGRAM ERWF

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum