text:
C    %W% %G%
C****************************************************************
C
C     File: ld_geldc.f
C
C     Purpose: Routine to load GE mulit-terminal d-c branch data from
C              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_ge.f
C
C****************************************************************
      integer function ld_geldc (xbuf, file, options, error)
      integer file, error, options(*)
      character *(*) xbuf

      include 'ipfinc/parametr.inc'

      include 'ipfinc/prt.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/alt_case.inc'
      include 'ipfinc/alpha.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/branch.inc'

      common /old_branch/ old_k1, old_k2, old_sect, old_type, old_id
      integer old_k1, old_k2, old_sect, old_type
      character old_id * 1

      character word(100)*10, id*1, getid*1
      integer find_bus, fnd_ptin, ftn_atoi, status, bptr1, bptr2,
     &        oldptr1, oldptr2, gtnwdcno, read_ge_file

      data oldptr1, oldptr2 / 0, 0/

      save oldptr1, oldptr2

      ld_geldc = 0
      error = 0

      last = lastch0 (xbuf)
      if (xbuf(last:last) .eq. '/') then
        status = read_ge_file (file, xbuf(last:))
        if (status .eq. 0) go to 340
        last = lastch0 (xbuf)
      endif
      call uscan (xbuf(1:last), word, nwrd, '~',  ' ')

      do i = nwrd+1, 45
        word(i) = ' '
      enddo

      num1 = ftn_atoi(word(1))
      num2 = ftn_atoi(word(4))

      id = getid (word(9))
      if (word(11) .eq. ' ') then
        status = 1
      else
        status = ftn_atoi(word(10))
      endif
      if (num1 .eq. 0) then
        if (ichar (xbuf(1:1)) .ge. ichar ('a') .and.
     &      ichar (xbuf(1:1)) .le. ichar ('z')) ld_geldc = 1
      else if (num2 .eq. 0 .or. status .eq. 0) then
      else
        if (ltot+1 .ge. MAXBRN) then
          write (errbuf(1), 10000) MAXBRN
10000     format ('More than ',i5,
     &            ' branch records. Overflow occurred at branch:')
          write (errbuf(2), 10010) xbuf(1:80)
10010     format(11x, '(', a, ')')
          call prterx ('E',2)
          ltot = 1
          error = 1
          go to 900
        endif
        npti1 = gtnwdcno (num1)
        npti2 = gtnwdcno (num2)
        kp1 = fnd_ptin (npti1)
        if (kp1 .le. 0) then
          write (errbuf(1), 10020) num1, num2
10020     format (' Bus1 on d-c branch record (', i6, 1x, i6,
     &      ') is not a valid bus number.')
          call prterx ('W',1)
          error = 1
          go to 900
        endif
        kp2 = fnd_ptin (npti2)
        if (kp2 .le. 0) then
          write (errbuf(1), 10030) num1, num2
10030     format (' Bus2 on d-c branch record (', i6, 1x, i6,
     &      ') is not a valid bus number.')
          call prterx ('W',1)
          error = 1
          go to 900
        endif
        k1 = find_bus (pti_name(kp1), pti_base(kp1))
        k2 = find_bus (pti_name(kp2), pti_base(kp2))
        if (k1 .le. 0) then
          write (errbuf(1), 10040) npti1, pti_name(kp1),
     &      pti_base(kp1)
10040     format (' D-C bus1 (', i6, 1x, a8, f7.1, ') is not in system.
     &')
          call prterx ('W',1)
          error = 1
          go to 900
        else if (k2 .le. 0) then
          write (errbuf(1), 10050) npti2, pti_name(kp2),
     &       pti_base(kp2)
10050     format (' D-C bus2 (', i6, 1x, a8, f7.1, ') is not in system.
     &')
          call prterx ('W',1)
          error = 1
          go to 900
        endif
        ltot = ltot + 1
        kbrnch(1,ltot) = 2
        do i = 2, 18
          kbrnch(i,ltot) = 0
        enddo
        call putchr (3, word(30), kbrnch(3,ltot))
        brnch(5,ltot) = ftn_atof(word(13))   ! Rdc
        brnch(6,ltot) = ftn_atof(word(14))   ! Ldc
        brnch(7,ltot) = ftn_atof(word(15))   ! Cdc

        const = 1000.0 * base(k1)
        ratea = const * ftn_atof(word(16)) ! Rate A
        rateb = const * ftn_atof(word(17)) ! Rate B
        ratec = const * ftn_atof(word(18)) ! Rate C
        rated = const * ftn_atof(word(19)) ! Rate D
        ratea = amin1 (ratea, 9999.0)
        rateb = amin1 (rateb, 9999.0)
        rateb = amin1 (rateb, 9999.0)
        if (options(5) .eq. 3) then
          brnch(4,ltot) = ratec            ! Rate
        else
          brnch(4,ltot) = ratea            ! Rate
        endif
C
C       Add extended ratings.
C
        if (options(6) .eq. 3) then
          rateln(1,ltot) = ratec            ! Rate C
        else
          rateln(1,ltot) = ratea            ! Rate A
        endif
        if (options(7) .eq. 1) then
          rateln(2,ltot) = ratea            ! Rate A
        else
          rateln(2,ltot) = ratec            ! Rate C
        endif

c       link up branch

        call lkbrdata (ltot, k1, k2, bptr1, bptr2, error)
        kx(bptr1) = k1

        ky(bptr1) = k2
        brid(bptr1) = id
        brsect(bptr1) = 0
        brtype(bptr1) = kbrnch(1,ltot)

        kx(bptr2) = k2
        ky(bptr2) = k1
        brid(bptr2) = id
        brsect(bptr2) = 0
        brtype(bptr2) = kbrnch(1,ltot)

        old_k1 = k1
        old_k2 = k2
        old_id = id
        old_sect = 0

      endif

      go to 900

  340 continue
      write (errbuf(1), 10060)
10060 format (' Premature E-O-F encountered reading Branch Data segment'
     &)
      write (errbuf(2), 10010) xbuf(1:80)
      call prterx ('W',2)
      ld_geldc = 1
      error = 1

  900 continue
      return
      end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
* From arpa!sol-michael.stanford.edu!mike 5 May 89 23:53:00 PDT
      SUBROUTINE LSQR  ( M, N, APROD, DAMP,
     $                   LENIW, LENRW, IW, RW,
     $                   U, V, W, X, SE,
     $                   ATOL, BTOL, CONLIM, ITNLIM, NOUT,
     $                   ISTOP, ITN, ANORM, ACOND, RNORM, ARNORM, XNORM)

      EXTERNAL           APROD
      INTEGER            M, N, LENIW, LENRW, ITNLIM, NOUT, ISTOP, ITN
      INTEGER            IW(LENIW)
      DOUBLE PRECISION   RW(LENRW), U(M), V(N), W(N), X(N), SE(N),
     $                   ATOL, BTOL, CONLIM, DAMP,
     $                   ANORM, ACOND, RNORM, ARNORM, XNORM
*-----------------------------------------------------------------------
*
*     LSQR  finds a solution x to the following problems:
*
*     1. Unsymmetric equations --    solve  A*x = b
*
*     2. Linear least squares  --    solve  A*x = b
*                                    in the least-squares sense
*
*     3. Damped least squares  --    solve  (   A    )*x = ( b )
*                                           ( damp*I )     ( 0 )
*                                    in the least-squares sense
*
*     where A is a matrix with m rows and n columns, b is an
*     m-vector, and damp is a scalar.  (All quantities are real.)
*     The matrix A is intended to be large and sparse.  It is accessed
*     by means of subroutine calls of the form
*
*                CALL APROD ( mode, m, n, x, y, LENIW, LENRW, IW, RW )
*
*     which must perform the following functions:
*
*                If MODE = 1, compute  y = y + A*x.
*                If MODE = 2, compute  x = x + A(transpose)*y.
*
*     The vectors x and y are input parameters in both cases.
*     If  mode = 1,  y should be altered without changing x.
*     If  mode = 2,  x should be altered without changing y.
*     The parameters LENIW, LENRW, IW, RW may be used for workspace
*     as described below.
*
*     The rhs vector b is input via U, and subsequently overwritten.
*
*
*     Note:  LSQR uses an iterative method to approximate the solution.
*     The number of iterations required to reach a certain accuracy
*     depends strongly on the scaling of the problem.  Poor scaling of
*     the rows or columns of A should therefore be avoided where
*     possible.
*
*     For example, in problem 1 the solution is unaltered by
*     row-scaling.  If a row of A is very small or large compared to
*     the other rows of A, the corresponding row of ( A  b ) should be
*     scaled up or down.
*
*     In problems 1 and 2, the solution x is easily recovered
*     following column-scaling.  Unless better information is known,
*     the nonzero columns of A should be scaled so that they all have
*     the same Euclidean norm (e.g., 1.0).
*
*     In problem 3, there is no freedom to re-scale if damp is
*     nonzero.  However, the value of damp should be assigned only
*     after attention has been paid to the scaling of A.
*
*     The parameter damp is intended to help regularize
*     ill-conditioned systems, by preventing the true solution from
*     being very large.  Another aid to regularization is provided by
*     the parameter ACOND, which may be used to terminate iterations
*     before the computed solution becomes very large.
*
*
*     Notation
*     --------
*
*     The following quantities are used in discussing the subroutine
*     parameters:
*
*     Abar   =  (   A    ),          bbar  =  ( b )
*               ( damp*I )                    ( 0 )
*
*     r      =  b  -  A*x,           rbar  =  bbar  -  Abar*x
*
*     rnorm  =  sqrt( norm(r)**2  +  damp**2 * norm(x)**2 )
*            =  norm( rbar )
*
*     RELPR  =  the relative precision of floating-point arithmetic
*               on the machine being used.  For example, on the IBM 370,
*               RELPR is about 1.0E-6 and 1.0D-16 in single and double
*               precision respectively.
*
*     LSQR  minimizes the function rnorm with respect to x.
*
*
*     Parameters
*     ----------
*
*     M       input      m, the number of rows in A.
*
*     N       input      n, the number of columns in A.
*
*     APROD   external   See above.
*
*     DAMP    input      The damping parameter for problem 3 above.
*                        (DAMP should be 0.0 for problems 1 and 2.)
*                        If the system A*x = b is incompatible, values
*                        of DAMP in the range 0 to sqrt(RELPR)*norm(A)
*                        will probably have a negligible effect.
*                        Larger values of DAMP will tend to decrease
*                        the norm of x and reduce the number of
*                        iterations required by LSQR.
*
*                        The work per iteration and the storage needed
*                        by LSQR are the same for all values of DAMP.
*
*     LENIW   input      The length of the workspace array IW.
*     LENRW   input      The length of the workspace array RW.
*     IW      workspace  An integer array of length LENIW.
*     RW      workspace  A real array of length LENRW.
*
*             Note:  LSQR  does not explicitly use the previous four
*             parameters, but passes them to subroutine APROD for
*             possible use as workspace.  If APROD does not need
*             IW or RW, the values LENIW = 1 or LENRW = 1 should
*             be used, and the actual parameters corresponding to
*             IW or RW  may be any convenient array of suitable type.
*
*     U(M)    input      The rhs vector b.  Beware that U is
*                        over-written by LSQR.
*
*     V(N)    workspace
*     W(N)    workspace
*
*     X(N)    output     Returns the computed solution x.
*
*     SE(N)   output     Returns standard error estimates for the
*                        components of X.  For each i, SE(i) is set
*                        to the value  rnorm * sqrt( sigma(i,i) / T ),
*                        where sigma(i,i) is an estimate of the i-th
*                        diagonal of the inverse of Abar(transpose)*Abar
*                        and  T = 1      if  m .le. n,
*                             T = m - n  if  m .gt. n  and  damp = 0,
*                             T = m      if  damp .ne. 0.
*
*     ATOL    input      An estimate of the relative error in the data
*                        defining the matrix A.  For example,
*                        if A is accurate to about 6 digits, set
*                        ATOL = 1.0E-6 .
*
*     BTOL    input      An extimate of the relative error in the data
*                        defining the rhs vector b.  For example,
*                        if b is accurate to about 6 digits, set
*                        BTOL = 1.0E-6 .
*
*     CONLIM  input      An upper limit on cond(Abar), the apparent
*                        condition number of the matrix Abar.
*                        Iterations will be terminated if a computed
*                        estimate of cond(Abar) exceeds CONLIM.
*                        This is intended to prevent certain small or
*                        zero singular values of A or Abar from
*                        coming into effect and causing unwanted growth
*                        in the computed solution.
*
*                        CONLIM and DAMP may be used separately or
*                        together to regularize ill-conditioned systems.
*
*                        Normally, CONLIM should be in the range
*                        1000 to 1/RELPR.
*                        Suggested value:
*                        CONLIM = 1/(100*RELPR)  for compatible systems,
*                        CONLIM = 1/(10*sqrt(RELPR)) for least squares.
*
*             Note:  If the user is not concerned about the parameters
*             ATOL, BTOL and CONLIM, any or all of them may be set
*             to zero.  The effect will be the same as the values
*             RELPR, RELPR and 1/RELPR respectively.
*
*     ITNLIM  input      An upper limit on the number of iterations.
*                        Suggested value:
*                        ITNLIM = n/2   for well-conditioned systems
*                                       with clustered singular values,
*                        ITNLIM = 4*n   otherwise.
*
*     NOUT    input      File number for printed output.  If positive,
*                        a summary will be printed on file NOUT.
*
*     ISTOP   output     An integer giving the reason for termination:
*
*                0       x = 0  is the exact solution.
*                        No iterations were performed.
*
*                1       The equations A*x = b are probably
*                        compatible.  Norm(A*x - b) is sufficiently
*                        small, given the values of ATOL and BTOL.
*
*                2       The system A*x = b is probably not
*                        compatible.  A least-squares solution has
*                        been obtained that is sufficiently accurate,
*                        given the value of ATOL.
*
*                3       An estimate of cond(Abar) has exceeded
*                        CONLIM.  The system A*x = b appears to be
*                        ill-conditioned.  Otherwise, there could be an
*                        error in subroutine APROD.
*
*                4       The equations A*x = b are probably
*                        compatible.  Norm(A*x - b) is as small as
*                        seems reasonable on this machine.
*
*                5       The system A*x = b is probably not
*                        compatible.  A least-squares solution has
*                        been obtained that is as accurate as seems
*                        reasonable on this machine.
*
*                6       Cond(Abar) seems to be so large that there is
*                        no point in doing further iterations,
*                        given the precision of this machine.
*                        There could be an error in subroutine APROD.
*
*                7       The iteration limit ITNLIM was reached.
*
*     ITN     output     The number of iterations performed.
*
*     ANORM   output     An estimate of the Frobenius norm of  Abar.
*                        This is the square-root of the sum of squares
*                        of the elements of Abar.
*                        If DAMP is small and if the columns of A
*                        have all been scaled to have length 1.0,
*                        ANORM should increase to roughly sqrt(n).
*                        A radically different value for ANORM may
*                        indicate an error in subroutine APROD (there
*                        may be an inconsistency between modes 1 and 2).
*
*     ACOND   output     An estimate of cond(Abar), the condition
*                        number of Abar.  A very high value of ACOND
*                        may again indicate an error in APROD.
*
*     RNORM   output     An estimate of the final value of norm(rbar),
*                        the function being minimized (see notation
*                        above).  This will be small if A*x = b has
*                        a solution.
*
*     ARNORM  output     An estimate of the final value of
*                        norm( Abar(transpose)*rbar ), the norm of
*                        the residual for the usual normal equations.
*                        This should be small in all cases.  (ARNORM
*                        will often be smaller than the true value
*                        computed from the output vector X.)
*
*     XNORM   output     An estimate of the norm of the final
*                        solution vector X.
*
*
*     Subroutines and functions used
*     ------------------------------
*
*     USER               APROD
*     BLAS               DCOPY, DNRM2, DSCAL (see Lawson et al. below)
*
*
*     Precision
*     ---------
*
*     The number of iterations required by LSQR will usually decrease
*     if the computation is performed in higher precision.  To convert
*     LSQR between single and double precision, change the words
*                        DOUBLE PRECISION
*                        DCOPY, DNRM2, DSCAL
*     to the appropriate FORTRAN and BLAS equivalents.
*     Also change 'D+' or 'E+' in the PARAMETER statement.
*
*
*     References
*     ----------
*
*     C.C. Paige and M.A. Saunders,  LSQR: An algorithm for sparse
*          linear equations and sparse least squares,
*          ACM Transactions on Mathematical Software 8, 1 (March 1982),
*          pp. 43-71.
*
*     C.C. Paige and M.A. Saunders,  Algorithm 583, LSQR: Sparse
*          linear equations and least-squares problems,
*          ACM Transactions on Mathematical Software 8, 2 (June 1982),
*          pp. 195-209.
*
*     C.L. Lawson, R.J. Hanson, D.R. Kincaid and F.T. Krogh,
*          Basic linear algebra subprograms for Fortran usage,
*          ACM Transactions on Mathematical Software 5, 3 (Sept 1979),
*          pp. 308-323 and 324-325.
*-----------------------------------------------------------------------
*
*
*     LSQR development:
*     22 Feb 1982: LSQR sent to ACM TOMS to become Algorithm 583.
*     15 Sep 1985: Final F66 version.  LSQR sent to "misc" in netlib.
*     13 Oct 1987: Bug (Robert Davies, DSIR).  Have to delete
*                     IF ( (ONE + DABS(T)) .LE. ONE ) GO TO 200
*                  from loop 200.  The test was an attempt to reduce
*                  underflows, but caused W(I) not to be updated.
*     17 Mar 1989: First F77 version.
*     04 May 1989: Bug (David Gay, AT&T).  When the second BETA is zero,
*                  RNORM = 0 and
*                  TEST2 = ARNORM / (ANORM * RNORM) overflows.
*                  Fixed by testing for RNORM = 0.
*     05 May 1989: Sent to "misc" in netlib.
*
*     Michael A. Saunders            (na.saunders @ NA-net.stanford.edu)
*     Department of Operations Research
*     Stanford University
*     Stanford, CA 94305-4022.
*-----------------------------------------------------------------------

*     Intrinsics and local variables

      INTRINSIC          ABS, MOD, SQRT
      INTEGER            I, NCONV, NSTOP
      DOUBLE PRECISION   DNRM2
      DOUBLE PRECISION   ALFA, BBNORM, BETA, BNORM,
     $                   CS, CS1, CS2, CTOL, DAMPSQ, DDNORM, DELTA,
     $                   GAMMA, GAMBAR, PHI, PHIBAR, PSI,
     $                   RES1, RES2, RHO, RHOBAR, RHBAR1, RHBAR2,
     $                   RHS, RTOL, SN, SN1, SN2,
     $                   T, TAU, TEST1, TEST2, TEST3,
     $                   THETA, T1, T2, T3, XXNORM, Z, ZBAR

      DOUBLE PRECISION   ZERO,           ONE
      PARAMETER        ( ZERO = 0.0D+0,  ONE = 1.0D+0 )

C      CHARACTER*16       ENTER, EXIT
C      CHARACTER*60       MSG(0:7)

C      DATA               ENTER /' Enter LSQR.    '/,
C     $                   EXIT  /' Exit  LSQR.    '/

C      DATA               MSG
C     $ / 'The exact solution is  X = 0',
C     $   'Ax - b is small enough, given ATOL, BTOL',
C     $   'The least-squares solution is good enough, given ATOL',
C     $   'The estimate of cond(Abar) has exceeded CONLIM',
C     $   'Ax - b is small enough for this machine',
C     $   'The least-squares solution is good enough for this machine',
C     $   'Cond(Abar) seems to be too large for this machine',
C     $   'The iteration limit has been reached' /
*-----------------------------------------------------------------------

*     Initialize.

C      IF (NOUT .GT. 0)
C     $   WRITE(NOUT, 1000) ENTER, M, N, DAMP, ATOL, CONLIM, BTOL, ITNLIM
      ITN    =   0
      ISTOP  =   0
      NSTOP  =   0
      CTOL   =   ZERO
      IF (CONLIM .GT. ZERO) CTOL = ONE / CONLIM
      ANORM  =   ZERO
      ACOND  =   ZERO
      BBNORM =   ZERO
      DAMPSQ =   DAMP**2
      DDNORM =   ZERO
      RES2   =   ZERO
      XNORM  =   ZERO
      XXNORM =   ZERO
      CS2    = - ONE
      SN2    =   ZERO
      Z      =   ZERO

      DO 10  I = 1, N
         V(I)  =  ZERO
         X(I)  =  ZERO
        SE(I)  =  ZERO
   10 CONTINUE

*     Set up the first vectors U and V for the bidiagonalization.
*     These satisfy  BETA*U = b,  ALFA*V = A(transpose)*U.

      ALFA   =   ZERO
      BETA   =   DNRM2 ( M, U, 1 )

      IF (BETA .GT. ZERO) THEN
         CALL DSCAL ( M, (ONE / BETA), U, 1 )
         CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW )
         ALFA   =   DNRM2 ( N, V, 1 )
      END IF

      IF (ALFA .GT. ZERO) THEN
         CALL DSCAL ( N, (ONE / ALFA), V, 1 )
         CALL DCOPY ( N, V, 1, W, 1 )
      END IF

      ARNORM =   ALFA * BETA
      IF (ARNORM .EQ. ZERO) GO TO 800

      RHOBAR =   ALFA
      PHIBAR =   BETA
      BNORM  =   BETA
      RNORM  =   BETA

C      IF (NOUT   .GT.  0  ) THEN
C         IF (DAMPSQ .EQ. ZERO) THEN
C             WRITE(NOUT, 1200)
C         ELSE
C             WRITE(NOUT, 1300)
C         END IF
C         TEST1  = ONE
C         TEST2  = ALFA / BETA
C         WRITE(NOUT, 1500) ITN, X(1), RNORM, TEST1, TEST2
C         WRITE(NOUT, 1600)
C      END IF

*     ------------------------------------------------------------------
*     Main iteration loop.
*     ------------------------------------------------------------------
  100 ITN    = ITN + 1

*     Perform the next step of the bidiagonalization to obtain the
*     next  BETA, U, ALFA, V.  These satisfy the relations
*                BETA*U  =  A*V  -  ALFA*U,
*                ALFA*V  =  A(transpose)*U  -  BETA*V.

      CALL DSCAL ( M, (- ALFA), U, 1 )
      CALL APROD ( 1, M, N, V, U, LENIW, LENRW, IW, RW )
      BETA   =   DNRM2 ( M, U, 1 )
      BBNORM =   BBNORM  +  ALFA**2  +  BETA**2  +  DAMPSQ

      IF (BETA .GT. ZERO) THEN
         CALL DSCAL ( M, (ONE / BETA), U, 1 )
         CALL DSCAL ( N, (- BETA), V, 1 )
         CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW )
         ALFA   =   DNRM2 ( N, V, 1 )
         IF (ALFA .GT. ZERO) THEN
            CALL DSCAL ( N, (ONE / ALFA), V, 1 )
         END IF
      END IF

*     Use a plane rotation to eliminate the damping parameter.
*     This alters the diagonal (RHOBAR) of the lower-bidiagonal matrix.

      RHBAR2 = RHOBAR**2  +  DAMPSQ
      RHBAR1 = SQRT( RHBAR2 )
      CS1    = RHOBAR / RHBAR1
      SN1    = DAMP   / RHBAR1
      PSI    = SN1 * PHIBAR
      PHIBAR = CS1 * PHIBAR

*     Use a plane rotation to eliminate the subdiagonal element (BETA)
*     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.

      RHO    =   SQRT( RHBAR2  +  BETA**2 )
      CS     =   RHBAR1 / RHO
      SN     =   BETA   / RHO
      THETA  =   SN * ALFA
      RHOBAR = - CS * ALFA
      PHI    =   CS * PHIBAR
      PHIBAR =   SN * PHIBAR
      TAU    =   SN * PHI

*     Update  X, W  and the standard error estimates.

      T1     =   PHI   / RHO
      T2     = - THETA / RHO
      T3     =   ONE   / RHO

      DO 200  I =  1, N
         T      =  W(I)
         X(I)   =  T1*T  +  X(I)
         W(I)   =  T2*T  +  V(I)
         T      = (T3*T)**2
         SE(I)  =  T     +  SE(I)
         DDNORM =  T     +  DDNORM
  200 CONTINUE

*     Use a plane rotation on the right to eliminate the
*     super-diagonal element (THETA) of the upper-bidiagonal matrix.
*     Then use the result to estimate  norm(X).

      DELTA  =   SN2 * RHO
      GAMBAR = - CS2 * RHO
      RHS    =   PHI    - DELTA * Z
      ZBAR   =   RHS    / GAMBAR
      XNORM  =   SQRT( XXNORM    + ZBAR **2 )
      GAMMA  =   SQRT( GAMBAR**2 + THETA**2 )
      CS2    =   GAMBAR / GAMMA
      SN2    =   THETA  / GAMMA
      Z      =   RHS    / GAMMA
      XXNORM =   XXNORM + Z**2

*     Test for convergence.
*     First, estimate the norm and condition of the matrix  Abar,
*     and the norms of  rbar  and  Abar(transpose)*rbar.

      ANORM  =   SQRT( BBNORM )
      ACOND  =   ANORM * SQRT( DDNORM )
      RES1   =   PHIBAR**2
      RES2   =   RES2  +  PSI**2
      RNORM  =   SQRT( RES1 + RES2 )
      ARNORM =   ALFA  * ABS( TAU )

*     Now use these norms to estimate certain other quantities,
*     some of which will be small near a solution.

      TEST1  =   RNORM /  BNORM
      TEST2  =   ZERO
      IF (RNORM .GT. ZERO) TEST2 = ARNORM / (ANORM * RNORM)
      TEST3  =   ONE   /  ACOND
      T1     =   TEST1 / (ONE  +  ANORM * XNORM / BNORM)
      RTOL   =   BTOL  +  ATOL *  ANORM * XNORM / BNORM

*     The following tests guard against extremely small values of
*     ATOL, BTOL  or  CTOL.  (The user may have set any or all of
*     the parameters  ATOL, BTOL, CONLIM  to zero.)
*     The effect is equivalent to the normal tests using
*     ATOL = RELPR,  BTOL = RELPR,  CONLIM = 1/RELPR.

      T3     =   ONE + TEST3
      T2     =   ONE + TEST2
      T1     =   ONE + T1
      IF (ITN .GE. ITNLIM) ISTOP = 7
      IF (T3  .LE. ONE   ) ISTOP = 6
      IF (T2  .LE. ONE   ) ISTOP = 5
      IF (T1  .LE. ONE   ) ISTOP = 4

*     Allow for tolerances set by the user.

      IF (TEST3 .LE. CTOL) ISTOP = 3
      IF (TEST2 .LE. ATOL) ISTOP = 2
      IF (TEST1 .LE. RTOL) ISTOP = 1
*     ==================================================================

*     See if it is time to print something.

      IF (NOUT  .LE.  0       ) GO TO 600
      IF (N     .LE. 40       ) GO TO 400
      IF (ITN   .LE. 10       ) GO TO 400
      IF (ITN   .GE. ITNLIM-10) GO TO 400
      IF (MOD(ITN,10) .EQ. 0  ) GO TO 400
      IF (TEST3 .LE.  2.0*CTOL) GO TO 400
      IF (TEST2 .LE. 10.0*ATOL) GO TO 400
      IF (TEST1 .LE. 10.0*RTOL) GO TO 400
      IF (ISTOP .NE.  0       ) GO TO 400
      GO TO 600

*     Print a line for this iteration.

  400 IF (1 .EQ. 1) GO TO 600
C  400 WRITE(NOUT, 1500) ITN, X(1), RNORM, TEST1, TEST2, ANORM, ACOND
C      IF (MOD(ITN,10) .EQ. 0) WRITE(NOUT, 1600)
*     ==================================================================

*     Stop if appropriate.
*     The convergence criteria are required to be met on  NCONV
*     consecutive iterations, where  NCONV  is set below.
*     Suggested value:  NCONV = 1, 2  or  3.

  600 IF (ISTOP .EQ. 0) NSTOP = 0
      IF (ISTOP .EQ. 0) GO TO 100
      NCONV  =   1
      NSTOP  =   NSTOP + 1
      IF (NSTOP .LT. NCONV  .AND.  ITN .LT. ITNLIM) ISTOP = 0
      IF (ISTOP .EQ. 0) GO TO 100
*     ------------------------------------------------------------------
*     End of iteration loop.
*     ------------------------------------------------------------------

*     Finish off the standard error estimates.

      T    =   ONE
      IF (M      .GT.   N )  T = M - N
      IF (DAMPSQ .GT. ZERO)  T = M
      T    =   RNORM / SQRT( T )

      DO 700  I = 1, N
         SE(I)  = T * SQRT( SE(I) )
  700 CONTINUE

*     Print the stopping condition.

  800 IF (1 .EQ. 1) GO TO 900
C  800 IF (NOUT .GT. 0) THEN
C         WRITE(NOUT, 2000) EXIT, ISTOP, ITN,
C     $                     EXIT, ANORM, ACOND,
C     $                     EXIT, RNORM, ARNORM,
C     $                     EXIT, BNORM, XNORM
C         WRITE(NOUT, 3000) EXIT, MSG(ISTOP)
C      END IF

  900 RETURN

*     ------------------------------------------------------------------
 1000 FORMAT(// 1P, A, '  Least-squares solution of  A*x = b'
     $    / ' The matrix  A  has', I7, ' rows   and', I7, ' columns'
     $    / ' The damping parameter is         DAMP   =', E10.2
     $    / ' ATOL   =', E10.2, 15X,        'CONLIM =', E10.2
     $    / ' BTOL   =', E10.2, 15X,        'ITNLIM =', I10)
 1200 FORMAT(// '   Itn       x(1)           Function',
     $   '     Compatible   LS        Norm A    Cond A' /)
 1300 FORMAT(// '   Itn       x(1)           Function',
     $   '     Compatible   LS     Norm Abar Cond Abar' /)
 1500 FORMAT(1P, I6, 2E17.9, 4E10.2)
 1600 FORMAT(1X)
 2000 FORMAT(/ 1P, A, 6X, 'ISTOP =', I3,   16X, 'ITN    =', I9
     $       /     A, 6X, 'ANORM =', E13.5, 6X, 'ACOND  =', E13.5
     $       /     A, 6X, 'RNORM =', E13.5, 6X, 'ARNORM =', E13.5,
     $       /     A, 6X, 'BNORM =', E13.5, 6X, 'XNORM  =', E13.5)
 3000 FORMAT( A, 6X, A )
*     ------------------------------------------------------------------
*     End of LSQR
      END

max_stars_repo_name: jcfr/ITK

max_issues_repo_name: jcfr/ITK

max_forks_repo_name: jcfr/ITK


! -----------------------------------------------------------------------------
!
! Copyright (c) 2017 Sam Cox, Roberto Sommariva
!
! This file is part of the AtChem2 software package.
!
! This file is covered by the MIT license which can be found in the file
! LICENSE.md at the top level of the AtChem2 distribution.
!
! -----------------------------------------------------------------------------

module solar_test
  use fruit
  use types_mod
  use solar_functions_mod
  implicit none

contains

  subroutine test_calcTheta
    use date_mod, only : currentYear, currentDayOfYear
    real(kind=DP) :: theta, pi, threshold

    threshold = 1.0e-8_DP
    pi = 4.0_DP * atan( 1.0_DP )

    currentYear = 2000_DI
    currentDayOfYear = 0_DI
    theta = calcTheta()
    call assert_true( theta == 0.0_DP, "calcTheta(), first day of 2000" )
    currentDayOfYear = 365_DI
    theta = calcTheta()
    call assert_true( theta == 2.0_DP*pi*365_DI/366_DI, "calcTheta(), last day of 2000" )

    currentYear = 2001_DI
    currentDayOfYear = 0_DI
    theta = calcTheta()
    call assert_true( theta == 0.0_DP, "calcTheta(), first day of 2001" )
    currentDayOfYear = 364_DI
    theta = calcTheta()
    call assert_true( abs(theta - 2.0_DP*pi*364_DI/365_DI) <= threshold, "calcTheta(), last day of 2001" )

    currentYear = 2004_DI
    currentDayOfYear = 0_DI
    theta = calcTheta()
    call assert_true( theta == 0.0_DP, "calcTheta(), first day of 2004" )
    currentDayOfYear = 365_DI
    theta = calcTheta()
    call assert_true( theta == 2.0_DP*pi*365_DI/366_DI, "calcTheta(), last day of 2004" )
  end subroutine test_calcTheta

  subroutine test_decFromTheta
    real(kind=DP) :: dec, pi, theta

    pi = 4.0_DP * atan( 1.0_DP )

    theta = 0.0_DP
    dec = decFromTheta( theta )
    call assert_true( dec == -0.402449_DP, "decFromTheta( theta ), theta = 0" )

    theta = pi
    dec = decFromTheta( theta )
    call assert_true( dec == 0.40276899999999999_DP, "decFromTheta( theta ), theta = pi" )

    theta = pi/2.0_DP
    dec = decFromTheta( theta )
    call assert_true( dec == 8.2452999999999985E-002_DP, "decFromTheta( theta ), theta = pi/2" )

    theta = pi*3.0_DP/2.0_DP
    dec = decFromTheta( theta )
    call assert_true( dec == -5.5100999999999921E-002_DP, "decFromTheta( theta ), theta = 3pi/2" )
  end subroutine test_decFromTheta

  subroutine test_calcDec
    use types_mod
    use date_mod, only : currentYear, currentDayOfYear
    real(kind=DP) :: dec0, dec364over365, dec365over366, threshold

    threshold = 1.0e-6_DP
    dec0 = 0.006918_DP - 0.399912_DP - 0.006758_DP - 0.002697_DP
    dec364over365 = -0.40369912461219781
    dec365over366 = -0.40369589165596537

    currentYear = 2000_DI
    currentDayOfYear = 0_DI
    call assert_true( calcDec() == dec0, "calcDec(), first day of 2000" )
    currentDayOfYear = 365_DI
    call assert_true( abs( calcDec() - dec365over366 ) < threshold, "calcDec(), last day of 2000" )

    currentYear = 2001_DI
    currentDayOfYear = 0_DI
    call assert_true( calcDec() == dec0, "calcDec(), first day of 2001" )
    currentDayOfYear = 364_DI
    call assert_true( abs( calcDec() - dec364over365 ) < threshold, "calcDec(), last day of 2001" )

    currentYear = 2004_DI
    currentDayOfYear = 0_DI
    call assert_true( calcDec() == dec0, "calcDec(), first day of 2004" )
    currentDayOfYear = 365_DI
    call assert_true( abs( calcDec() - dec365over366 ) < threshold, "calcDec(), last day of 2004" )
  end subroutine test_calcDec

  subroutine test_calcEQT
    use types_mod
    use date_mod, only : currentYear, currentDayOfYear
    real(kind=DP) :: eqt0, threshold

    threshold = 1.0e-8_DP
    eqt0 = 0.000075_DP + 0.001868_DP - 0.014615_DP

    currentYear = 2000_DI
    currentDayOfYear = 0_DI
    call assert_true( abs( calcEQT() - eqt0 ) < threshold, "calcEQT(), first day of 2000" )
    currentDayOfYear = 365_DI
    call assert_true( abs( calcEQT() + 1.0710769160677822E-002 ) < threshold, "calcEQT(), last day of 2000" )

    currentYear = 2001_DI
    currentDayOfYear = 0_DI
    call assert_true( abs( calcEQT() - eqt0 ) < threshold, "calcEQT(), first day of 2001" )
    currentDayOfYear = 364_DI
    call assert_true( abs( calcEQT() + 1.0705374687579837E-002 ) < threshold, "calcEQT(), last day of 2001" )

    currentYear = 2004_DI
    currentDayOfYear = 0_DI
    call assert_true( abs( calcEQT() - eqt0 ) < threshold, "calcEQT(), first day of 2004" )
    currentDayOfYear = 365_DI
    call assert_true( abs( calcEQT() + 1.0710769160677822E-002 ) < threshold, "calcEQT(), last day of 2004" )
  end subroutine test_calcEQT

  ! TODO:  subroutine test_calcZenith

end module solar_test


text:
Peter J. Shields Avenue, or Shields Ave. for short, is the default address for the university, One Shields Avenue. The street runs west from the corner of 1st Street 1st and A Street A until it Ts into California Avenue. It passes the north side of Shields Library.  Try as you like, you will not find an actual, 1 Shields Avenue address to go with a building.

The Avenue is named after Peter J. Shields, a man who helped create UC Davis.

Other things on Shields Ave.:
  Olson Hall
  Walker Hall
  AOB IV
  Cross Cultural Center
  Hart Hall

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


          program demo_rd
          use M_io, only : rd
          implicit none
          character(len=:),allocatable :: mystring
          doubleprecision              :: d
          real                         :: r
          integer                      :: i
          logical                      :: l

          INFINITE: do
             mystring=rd('Enter string or "STOP":',default='Today')
             if(mystring.eq.'STOP')stop
             i=rd('Enter integer:',default=huge(0))
             r=rd('Enter real:',default=huge(0.0))
             d=rd('Enter double:',default=huge(0.0d0))
             l=rd('Enter logical:',default=.false.)

             write(*,*)'I=', i, 'R=', r, 'D=',d,  'MYSTRING=', mystring
             write(*,*)'L=', l
          enddo INFINITE

          end program demo_rd


      subroutine nut_nlch

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine simulates the loss of nitrate via surface runoff,
!!    lateral flow, tile flow, and percolation out of the profile

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    nperco      |none          |nitrate percolation coefficient (0-1)
!!                               |0:concentration of nitrate in surface runoff
!!                               |  is zero
!!                               |1:surface runoff has same concentration of
!!                               |  nitrate as percolate
!!    surfq(:)    |mm H2O        |surface runoff generated on day in HRU
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Exp, Max, Min

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use basin_module
      use organic_mineral_mass_module
      use hru_module, only : hru, latno3, percn, surqno3, tileno3, surfq, ihru, qtile
      use soil_module

      implicit none

      integer :: j         !none          |HRU number
      integer :: jj        !none          |counter
      real :: sro          !mm H2O        |surface runoff
      real :: ssfnlyr      !kg N/ha       |nitrate transported in lateral flow from layer
      real :: percnlyr     !kg N/ha       |nitrate leached to next lower layer with
                           !              |percolation
      real :: vv           !mm H2O        |water mixing with nutrient in layer
      real :: vno3         !              |
      real :: co           !kg N/mm       |concentration of nitrate in solution
      real :: cosurf       !kg N/mm       |concentration of nitrate in surface runoff
      real :: nloss        !frac          |nloss based on half life
      real :: ww           !varies        |variable to hold intermediate calculation
                           !              |result

      j = ihru

      percnlyr = 0.

      do jj = 1, soil(j)%nly

        !! add nitrate leached from layer above
        soil1(j)%mn(jj)%no3 = soil1(j)%mn(jj)%no3 + percnlyr
	  if (soil1(j)%mn(jj)%no3 < 1.e-6) soil1(j)%mn(jj)%no3 = 0.0

        !! determine concentration of nitrate in mobile water
        if (jj == 1) then
          sro = surfq(j)
        else
          sro = 0.
        end if
        vv = soil(j)%ly(jj)%prk + sro + soil(j)%ly(jj)%flat + 1.e-10
        if (hru(j)%lumv%ldrain == jj) vv = vv + qtile
        ww = -vv / ((1. - soil(j)%anion_excl) * soil(j)%phys(jj)%ul)
        vno3 = soil1(j)%mn(jj)%no3 * (1. - Exp(ww))
        co = Max(vno3 / vv, 0.)

        !! calculate nitrate in surface runoff
        cosurf = bsn_prm%nperco * co
        if (jj == 1) then
          surqno3(j) = surfq(j) * cosurf
          surqno3(j) = Min(surqno3(j), soil1(j)%mn(jj)%no3)
          soil1(j)%mn(jj)%no3 = soil1(j)%mn(jj)%no3 - surqno3(j)
        endif
        !Daniel 1/2012
        !! calculate nitrate in tile flow
        if (hru(j)%lumv%ldrain == jj) then
           ! tileno3(j) = bsn_prm%nperco * co * qtile     !Daniel 1/2012
           tileno3(j) = co * qtile     !Daniel 1/2012
          tileno3(j) = Min(tileno3(j), soil1(j)%mn(jj)%no3)
          soil1(j)%mn(jj)%no3 = soil1(j)%mn(jj)%no3 - tileno3(j)
        end if
        !Daniel 1/2012

        !! calculate nitrate in lateral flow
        if (jj == 1) then
          ssfnlyr = cosurf * soil(j)%ly(jj)%flat
        else
          ssfnlyr = co * soil(j)%ly(jj)%flat
        end if
        ssfnlyr = Min(ssfnlyr, soil1(j)%mn(jj)%no3)
        latno3(j) = latno3(j) + ssfnlyr
        soil1(j)%mn(jj)%no3 = soil1(j)%mn(jj)%no3 - ssfnlyr

        !! calculate nitrate in percolate
        percnlyr = co * soil(j)%ly(jj)%prk
        percnlyr = Min(percnlyr, soil1(j)%mn(jj)%no3)
        soil1(j)%mn(jj)%no3 = soil1(j)%mn(jj)%no3 - percnlyr
      end do

      !! calculate nitrate leaching from soil profile
      percn(j) = percnlyr

      nloss = (2.18 * hru(j)%topo%dis_stream - 8.63) / 100.
      nloss = Max(0.,nloss)
      nloss = Amin1(1.,nloss)
      latno3(j) = (1. - nloss) * latno3(j)

      return
      end subroutine nut_nlch


text:
  subroutine binned_cross_section_Raman( &
    maxlay,           temp,             jmax,            &
    nsmax,            ns_elas,          wavelength_ns,   &
    wavenr_ns,        wavenrshift_ns,   dnu,             &
    relN2,            relO2,            relAr,           &
    csray_N2,         csray_O2,         csray_Ar,        &
    eps_N2,           eps_O2,           linestr_Sto_N2,  &
    linestr_aSt_N2,   linestr_Sto_O2,   linestr_aSt_O2,  &
    dwn_Sto_N2,       dwn_aSt_N2,       dwn_Sto_O2,      &
    dwn_aSt_O2,       csram_Sto_N2,     csram_aSt_N2,    &
    csram_Sto_O2,     csram_aSt_O2,     csram_N2,        &
    csram_O2,         csram_air,        scat_type,       &
    weightCabRam_air)

!-----------------------------------------------------------
! Subroutine:
!    csrambin
!
! Requires:
!    binner
!    shift_column
!
!
! Purpose:
!
! Count count_ns only for ns2=ns2check and k=kcheck
!------------------------------------------------------------------

  implicit none

  integer,parameter::dp=selected_real_kind(8) !double precision

  integer,intent(in):: &
    jmax,   &  !maximum rotational qunatum number J
    maxlay, &  !nr of layers
    nsmax,  &     !nr of wavelengths
    ns_elas

  real(kind=dp),intent(in):: &
    eps_N2(0:nsmax),   &
    eps_O2(0:nsmax),   &
    csray_N2(0:nsmax), &
    csray_O2(0:nsmax), &
    csray_Ar(0:nsmax)

! wavenrshifts
  real(kind=dp),intent(in):: &
    dwn_Sto_N2(0:jmax), &
    dwn_aSt_N2(0:jmax), &
    dwn_Sto_O2(0:jmax), &
    dwn_aSt_O2(0:jmax)

  real(kind=dp),dimension(maxlay,0:jmax,0:nsmax),intent(in):: &
    csram_Sto_N2, &
    csram_aSt_N2, &
    csram_Sto_O2, &
    csram_aSt_O2

! wavelengths around incident light (ns2_elas) [nm]
  real(kind=dp),intent(in):: &
    wavelength_ns(0:nsmax),  &
    wavenr_ns(0:nsmax),      &
    wavenrshift_ns(0:nsmax), &
    dnu

  real(kind=dp),intent(in):: &
     temp(maxlay), & !temperature [K] for each layer
     relN2, & !N2 volume mixing ratio
     relO2, & !O2 volume mixing ratio
     relAr

  real(kind=dp),dimension(0:jmax,maxlay),intent(in):: &
    linestr_Sto_N2, &
    linestr_aSt_N2, &
    linestr_Sto_O2, &
    linestr_aSt_O2

  integer::ns1,ns2,k,j,ns

  real(kind=dp):: &
    testfac

  integer:: &
    kcheck,  &
    ns2check, &
    ns1r,     &
    count_ns(0:nsmax)

  real(kind=dp),dimension(maxlay,0:nsmax,0:nsmax),intent(out):: &
    csram_N2, &
    csram_O2, &
    csram_air    !Raman cross section for air

  integer,intent(out):: &
    scat_type(0:nsmax,0:nsmax) !0=nothing; 1=Cab; 2=Ram; 3=Cab+ram

  real(kind=dp),intent(out):: &
    weightCabRam_air(maxlay,0:nsmax) !how many inel. light
    !added top Cabbannes scattering for each ns2

  logical,parameter:: &
    verbose = .false.,            &
    energy_conservation = .true.

!-----------------------------------------------------------------------
! Reset Raman cross section
!----------------------------------------------------------------------

  do ns2=0,nsmax
     do ns1=0,nsmax
        do k=1,maxlay
           csram_N2(k,ns1,ns2)  = 0.d0
           csram_O2(k,ns1,ns2)  = 0.d0
           csram_air(k,ns1,ns2) = 0.d0
        enddo
        scat_type(ns1,ns2) = 0
     enddo
  enddo

!-----------------------------------------------------------------------
! Bin Raman lines and count (non-zero) Raman lines per bin
!
! incoming wavelength ns_elas
!-----------------------------------------------------------------------

! select layer to count Raman lines;
! count_ns is determined for this layer
  kcheck = maxlay

! count_ns is determined for wavenumber index ns2check
  ns2check = 0

! reset counter
  do ns1=0,nsmax
     count_ns(ns1) = 0
  enddo

! loop over incident wavenr
  do ns2=0,nsmax

!    reset column
     do ns1=0,nsmax
        do k=1,maxlay
           csram_N2(k,ns1,ns2)  = 0.d0
           csram_O2(k,ns1,ns2)  = 0.d0
           csram_air(k,ns1,ns2) = 0.d0
         enddo
     enddo

!-------------------------------------------------------------
!    bin Raman lines per column
!-------------------------------------------------------------

!    Scattering type
     scat_type(0:nsmax,ns2) = 2 !Raman scattering
     scat_type(ns_elas,ns2) = 1 !Cabannes scattering

!    Initial value, pure Cabannes scattering for elements
!    (ns_elas,ns2)
     weightCabRam_air(1:maxlay,ns2) = 1.d0

     call bin_Raman_lines_per_column(  &
          ns2,            ns2check,      kcheck,         &
          jmax,           maxlay,        nsmax,          &
          ns_elas,        wavelength_ns, wavenr_ns,      &
          wavenrshift_ns, dnu,           relN2,          &
          relO2,          dwn_Sto_N2,    dwn_aSt_N2,     &
          dwn_Sto_O2,     dwn_aSt_O2,    csram_Sto_N2,   &
          csram_aSt_N2,   csram_Sto_O2,  csram_aSt_O2,   &
          csram_N2,       csram_O2,      csram_air,      &
          count_ns)

!--------------------------------------------------------------
!    Special test case
!--------------------------------------------------------------

     goto 444

     if(ns2.eq.ns_elas)then
        write(*,*) 'csram = 0.0 0.025 0.05 0.800 0.05 0.025 0.00'
     endif

     do k=1,maxlay

        csram_N2(k,0:nsmax,ns2)  = 0.d0
        csram_O2(k,0:nsmax,ns2)  = 0.d0
        csram_air(k,0:nsmax,ns2) = 0.d0

        do ns=-2,2
           if(abs(ns).eq.0) testfac = 0.800d0
           if(abs(ns).eq.1) testfac = 0.050d0
           if(abs(ns).eq.2) testfac = 0.025d0
           csram_N2(k,ns_elas+ns,ns2) = testfac*csray_N2(ns2)
           csram_O2(k,ns_elas+ns,ns2) = testfac*csray_O2(ns2)
           csram_air(k,ns_elas+ns,ns2) = &
              relN2*csram_N2(k,ns_elas+ns,ns2) + &
                relO2*csram_O2(k,ns_elas+ns,ns2)
        enddo

     enddo

 444 continue

     if(verbose)then
        do ns1=0,nsmax
           write(*,'(F8.2,I5,3D16.8,I4,f8.3)') &
             wavenrshift_ns(ns1),         &
             count_ns(ns1),               &
             csram_N2(kcheck,ns1,ns2),    &
             csram_O2(kcheck,ns1,ns2),    &
             csram_air(kcheck,ns1,ns2),   &
             scat_type(ns1,ns2),          &
             weightCabRam_air(kcheck,ns2)
        enddo
        write(*,*) '---'
      endif

!---------------------------------------------------------------
!    Energy conservation
!    Make sure that sum(csram)+cscab = csray
!    (might be obsolete)
!---------------------------------------------------------------

     if(energy_conservation)then

        call energy_conservation_column( &
          ns2,       ns_elas,   nsmax,        &
          maxlay,    relN2,     relO2,        &
          eps_N2,    eps_O2,    csray_N2,     &
          csray_O2,  csram_N2,  csram_O2,     &
          csram_air, scat_type,  weightCabRam_air)

      else

         write(*,*) 'csrambin.f: skip energy conservation'

      endif

!-------------------------------------------------------------
!    shift elements in column so that
!        1 1 1 1 1           3 2 1 0 0
!        2 2 2 2 2           4 3 2 1 0
!        3 3 3 3 3  becomes  5 4 3 2 1
!        4 4 4 4 4           0 5 4 3 2
!        5 5 5 5 5           0 0 5 4 3
!    and determine cabram weights (correct for border-effects)
!-------------------------------------------------------------

     call shift_elements_for_each_column( &
          ns2,       ns_elas,      nsmax,        &
          maxlay,    relN2,        relO2,        &
          relAr,     eps_N2,       eps_O2,       &
          csray_N2,  csray_O2,     csray_Ar,     &
          csram_N2,  csram_O2,     csram_air,    &
          scat_type,  weightCabRam_air)

     if(verbose)then
        do ns1=0,nsmax
           write(*,'(F8.2,I5,3D16.8,I4,F8.3)') &
             wavenrshift_ns(ns1),         &
             count_ns(ns1),               &
             csram_N2(kcheck,ns1,ns2),    &
             csram_O2(kcheck,ns1,ns2),    &
             csram_air(kcheck,ns1,ns2),   &
             scat_type(ns1,ns2),          &
             weightCabRam_air(kcheck,ns2)
        enddo
        write(*,*) '===================================---'
    endif

  enddo

  return
  end subroutine binned_cross_section_Raman

max_stars_repo_name: rutgervandeelen/das

max_issues_repo_name: rutgervandeelen/das

max_forks_repo_name: rutgervandeelen/das


text:
*
* $Id$
*

*    ************************************
*    *					*
*    *	    gen_PW92_unrestricted	*
*    *					*
*    ************************************
*
*    This function returns the PW92 LDA correlation
*  energy density, xce, and its derivatives with respect
*  to nup, ndn.
*
*   Entry - n2ft3d     : number of grid points
*           dn_in(*,2) : spin densites nup and ndn
*
*   Exit - xce(*)  : PBE96 energy density
*        - fn(*,2) : d(n*xce)/dnup, d(n*xce)/dndn

      subroutine gen_PW92_BW_unrestricted(n2ft3d,
     >                           dn_in,
     >                           xce,fn)
      implicit none

      integer n2ft3d
      real*8 dn_in(n2ft3d,2)
      real*8 xce(n2ft3d)
      real*8 fn(n2ft3d,2)

*     **** Density cutoff parameter ****
      real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2
      parameter (DNS_CUT = 1.0d-20)
      parameter (ETA=1.0d-20)
      parameter (ETA2=1.0d-14)
      parameter (alpha_zeta=(1.0d0-ETA2))
      parameter (alpha_zeta2=(1.0d0-ETA2))

c     ****** Perdew-Wang92 LDA correlation coefficients *******
      real*8 GAM,iGAM,FZZ,iFZZ
      parameter (GAM  	= 0.519842099789746329d0)
      parameter (iGAM  	= 1.0d0/GAM)
      parameter (FZZ    = (8.0d0/(9.0d0*GAM)) )
      parameter (iFZZ    = 0.125d0*9.0d0*GAM)

      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
      parameter (A_1  = 0.0310907d0)
      !parameter (A_1  = 0.031091d0)
      parameter (A1_1 =	0.2137000d0)
      parameter (B1_1 =	7.5957000d0)
      parameter (B2_1 =	3.5876000d0)
      parameter (B3_1 =	1.6382000d0)
      parameter (B4_1 =	0.4929400d0)

      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
      parameter (A_2  =  0.01554535d0)
      !parameter (A_2  =  0.015545d0)
      parameter (A1_2 =	 0.20548000d0)
      parameter (B1_2 =	14.11890000d0)
      parameter (B2_2 =	 6.19770000d0)
      parameter (B3_2 =	 3.36620000d0)
      parameter (B4_2 =	 0.62517000d0)

      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
      parameter (A_3  =  0.0168869d0)
      !parameter (A_3  =  0.016887d0)
      parameter (A1_3 =	 0.1112500d0)
      parameter (B1_3 =	10.3570000d0)
      parameter (B2_3 =	 3.6231000d0)
      parameter (B3_3 =	 0.8802600d0)
      parameter (B4_3 =	 0.4967100d0)

c     **** other constants ****
      real*8 onethird,fourthird
      real*8 onethirdm
      parameter (onethird=1.0d0/3.0d0)
      parameter (onethirdm=-1.0d0/3.0d0)
      parameter (fourthird=4.0d0/3.0d0)

c     **** local variables ****
      integer i
      real*8 n
      real*8 nup
      real*8 ndn
      real*8 n_onethird,pi,rs_scale
      real*8 rs         ! Wigner radius
      real*8 rss        ! rss  = sqrt(rs)
      real*8 rs_n       ! rs_n = n*drs/dn
      real*8 zet
      real*8 zet_nup    ! zet_nup = n*dzet/dnup
      real*8 zet_ndn    ! zet_nup = n*dzet/dnup
      real*8 zetp_1_3,zetm_1_3
      real*8 zetpm_1_3,zetmm_1_3
      real*8 FZ,z4
      real*8 F

      real*8 ecu,ecp,eca,ec,ec_lda
      real*8 ecu_rs,ecp_rs,eca_rs
      real*8 ec_lda_rs,ec_lda_zet  ! d(ec_lda)/drs, d(ec_lda)/dzet
      real*8 ec_lda_nup,ec_lda_ndn ! n*d(ec_lda)/dnup, n*d(ec_lda)/dndn
      real*8 fncup,fncdn           ! d(n*ec)/dnup, d(n*ec)/dndn

      pi = 4.0d0*datan(1.0d0)
      rs_scale = (0.75d0/pi)**onethird

!$OMP DO
      do i=1,n2ft3d
         nup     = dn_in(i,1)+ETA
         ndn     = dn_in(i,2)+ETA

c        *******************************************************************
c        ***** calculate polarized correlation energies and potentials *****
c        *******************************************************************
         n = nup+ndn
         zet = (nup-ndn)/n
         zet_nup = -(zet - 1.0d0)
         zet_ndn = -(zet + 1.0d0)
         zetpm_1_3 = (1.0d0+zet*alpha_zeta)**onethirdm
         zetmm_1_3 = (1.0d0-zet*alpha_zeta)**onethirdm
         zetp_1_3  = (1.0d0+zet*alpha_zeta)*zetpm_1_3**2
         zetm_1_3  = (1.0d0-zet*alpha_zeta)*zetmm_1_3**2

         F =(  (1.0d0+zet*alpha_zeta)*zetp_1_3
     >       + (1.0d0-zet*alpha_zeta)*zetm_1_3
     >       - 2.0d0)*iGAM

         FZ = (zetp_1_3 - zetm_1_3)*(alpha_zeta*fourthird*iGAM)

*        **** calculate Wigner radius ****
         rs    = rs_scale/(n**onethird)
         rss   = dsqrt(rs)

*        **** calculate n*drs/dn ****
c        rs_n = onethirdm*rs/n
         rs_n = onethirdm*rs

c        **************************************************
c        ***** compute LSDA correlation energy density ****
c        **************************************************
         call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ecu,ecu_rs)
         call LSDT(A_2,A1_2,B1_2,B2_2,B3_2,B4_2,rss,ecp,ecp_rs)
         call LSDT(A_3,A1_3,B1_3,B2_3,B3_3,B4_3,rss,eca,eca_rs)

         z4 = zet**4

         ec_lda = ecu*(1.0d0-F*z4)
     >          + ecp*F*z4
     >          - eca*F*(1.0d0-z4)/FZZ

         ec_lda_rs = ecu_rs*(1.0d0-F*z4)
     >             + ecp_rs*F*z4
     >             - eca_rs*F*(1.0d0-z4)/FZZ

         ec_lda_zet = (4.0d0*(zet**3)*F + FZ*z4)*(ecp-ecu+eca*iFZZ)
     >              - FZ*eca*iFZZ

         ec_lda_nup = ec_lda_zet
     >              - zet * ec_lda_zet
     >              + rs_n * ec_lda_rs
         ec_lda_ndn = -ec_lda_zet
     >              - zet  * ec_lda_zet
     >              + rs_n * ec_lda_rs

         ec = ec_lda

         fncup  = ec + (ec_lda_nup)
         fncdn  = ec + (ec_lda_ndn)

         xce(i)   = ec
         fn(i,1)  = fncup
         fn(i,2)  = fncdn
      end do
!$OMP END DO

      return
      end

*    ************************************
*    *					*
*    *	    gen_PW92_BW_restricted	*
*    *					*
*    ************************************
*
*   This routine calculates the PW92 LDA correlation
*   potential(xcp) and energy density(xce).
*
*
*   Entry - n2ft3d     : number of grid points
*           rho_in(*) :  density (nup+ndn)
*
*     Exit  - xce(n2ft3d) : PW92 correlation energy density
*             fn(n2ft3d)  : d(n*xce)/dn
*
      subroutine gen_PW92_BW_restricted(n2ft3d,rho_in,
     >                                xce,fn)
      implicit none

      integer    n2ft3d
      real*8     rho_in(n2ft3d)
      real*8     xce(n2ft3d)
      real*8     fn(n2ft3d)

*     **** Density cutoff parameter ****
      real*8 DNS_CUT,ETA
      parameter (DNS_CUT = 1.0d-20)
      parameter (ETA     = 1.0d-20)

c     ****** Perdew-Wang92 LDA correlation coefficients *******
      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
      parameter (A_1  = 0.0310907d0)
      parameter (A1_1 =	0.2137000d0)
      parameter (B1_1 =	7.5957000d0)
      parameter (B2_1 =	3.5876000d0)
      parameter (B3_1 =	1.6382000d0)
      parameter (B4_1 =	0.4929400d0)

      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
      parameter (A_2  =  0.01554535d0)
      parameter (A1_2 =	 0.20548000d0)
      parameter (B1_2 =	14.11890000d0)
      parameter (B2_2 =	 6.19770000d0)
      parameter (B3_2 =	 3.36620000d0)
      parameter (B4_2 =	 0.62517000d0)

      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
      parameter (A_3  =  0.0168869d0)
      parameter (A1_3 =	 0.1112500d0)
      parameter (B1_3 =	10.3570000d0)
      parameter (B2_3 =	 3.6231000d0)
      parameter (B3_3 =	 0.8802600d0)
      parameter (B4_3 =	 0.4967100d0)

c     **** other constants ****
      real*8 onethird
      parameter (onethird=1.0d0/3.0d0)

c     **** local variables ****
      integer i
      real*8 n
      real*8 n_onethird,pi,rs_scale
      real*8 rs,rss

      real*8 ec_lda
      real*8 ec_lda_rs
      real*8 ec
      real*8 fnc

      pi         = 4.0d0*datan(1.0d0)
      rs_scale   = (0.75d0/pi)**onethird

!$OMP DO
      do i=1,n2ft3d
         n     = rho_in(i)+ETA

*        *********************************************************************
c        ***** calculate unpolarized correlation energies and potentials *****
*        *********************************************************************

c        **** calculate rs and t ****
         rs    = rs_scale/(n**onethird)
         rss   = dsqrt(rs)

c        **** unpolarized LDA correlation energy ****
c        **** ec_p = correlation energy          ****
c        ****   ec_p_rs = dec_p/drs              ****
c        ****   uc_p    = dec_p/dn               ****
         call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ec_lda,ec_lda_rs)

         ec   = ec_lda
         fnc = ec  - (onethird*rs*ec_lda_rs)

         xce(i) = ec
         fn(i)  = fnc

      end do
!$OMP END DO

      return
      end

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
      SUBROUTINE NEWTON(TD4,TD5,TD6,IP,DTU)
*
*
*       Newton-Raphson iteration.
*       -------------------------
*
      INCLUDE 'common6.h'
      COMMON/SLOW0/  RANGE,ISLOW(10)
      PARAMETER (ONE30=1.0/30.0D0,ONE120=1.0/120.0D0,ZZ=1.0/720.0D0)
*
*
*       Include possible slow-down factor.
      IMOD = KSLOW(IP)
*     IF (IMOD.GT.1) THEN
          ZMOD = FLOAT(ISLOW(IMOD))
*     END IF
*
*       Include improved tolerance near small pericentre (R/a < 1.0D-02).
      IF (R(IP)*ABS(H(IP))/BODY(N+IP).LT.5.0D-03) THEN
          TOL = 1.0D-10
      ELSE
          TOL = 1.0D-08
      END IF
*
      DTU0 = DTAU(IP)
*       Note TD4 is 1/2 fourth derivative and TD5 is 1/4 fifth derivative.
      DO 10 I = 1,10
          FX = (((((ZZ*TD6*DTU + ONE30*TD5)*DTU + ONE12*TD4)*DTU +
     &             ONE6*TDOT3(IP))*DTU + 0.5*TDOT2(IP))*DTU + R(IP))*DTU
          FD = ((((ONE120*TD6*DTU + ONE6*TD5)*DTU + ONE3*TD4)*DTU +
     &             0.5*TDOT3(IP))*DTU + TDOT2(IP))*DTU + R(IP)
          IF (IMOD.GT.1) THEN
              FX = ZMOD*FX
              FD = ZMOD*FD
          END IF
          FX = FX - STEP(2*IP-1)
          DTU = DTU - FX/FD
          IF (DTU.LT.0.0D0) THEN
              DTU = 0.0
          ELSE IF (DTU.GT.DTU0) THEN
              DTU = DTU0
          END IF
*      Note last value of DTU ensures the actual tolerance is < 1D-12.
          IF (ABS(FX).LT.TOL.AND.DTU.GT.0.0D0) GO TO 30
   10 CONTINUE
*
      WRITE (6,20)  NSTEPU, FX, FD, DTU, DTU0
   20 FORMAT (' DANGER NEWTON!    # FX FD DTU DTU0 ',I11,1P,4E10.2)
      STOP
*
   30 RETURN
*
      END


text:
	module comflow
!     module comflow
!***********************************************************************
!  Copyright, 1993, 2004,  The  Regents of the University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
!D1
!D1  PURPOSE
!D1
!D1  Include file for array variables and pointers for particle tracking flow.
!D1
!***********************************************************************
!D2
!D2  REVISION HISTORY
!D2
!D2 $Log:   /pvcs.config/fehm90/src/comflow.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:42:30   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 08:56:48   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:05:34   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:22:32   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 11:56:52   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:39:42 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
!D2
!D2    Rev 1.4   Wed May 29 14:28:28 1996   hend
!D2 Added variable diffusion with water content
!D2
!D2    Rev 1.3   Thu Mar 21 13:20:08 1996   hend
!D2 Added variable for trac index ind. of istrw
!D2
!D2    Rev 1.2   Mon Mar 04 15:58:54 1996   hend
!D2 Removed uneccessary calculations from coneq1 and added trac input option
!D2
!D2    Rev 1.1   Tue Jan 16 15:55:14 1996   zvd
!D2 Added prolog.
!D2
!***********************************************************************
!D3
!D3  REQUIREMENTS TRACEABILITY
!D3
!D3  N/A
!D3
!***********************************************************************
!D4
!D4  SPECIAL COMMENTS AND REFERENCES
!D4
!D4  None
!D4
!***********************************************************************

	real*8, allocatable :: a_axy(:),a_vxy(:),a_wvxy(:),alphaconl(:)
	real*8, allocatable :: alphaconv(:),sehvell(:),sehvelv(:)
	real*8, allocatable :: sehdiff(:),sehdiffv(:)
!       real*8, allocatable :: thetaresid(:,:), thetaresidv(:,:)
        integer dispsame,sehdonevel,hvliquid,hvvapor,sehsize,numflux

c s kelkar 3 July 2014, for calculating heat flow vectors
	logical flag_heat_out
	real*8, allocatable :: e_axy_adv(:), e_axy_cond(:)
	real*8, allocatable :: e_adv_nodal(:,:),e_cond_nodal(:,:)

	end module comflow

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
!>-----------------------------------------------------------------------------------------------
!> @brief       Averaging techniques to obtain mean orbital elements
!!
!> @author      Vitali Braun (VB)
!> @author      Christopher Kebschull (CHK)
!!
!> @date

!!
VB: 30.07.2013 (initial design)

!!
VB: 16.01.2014 (now using 'units' module)

!!
CHK: 13.11.2015 (updated to use libslam, added some more documentation)

!!

!!
!> @details     This module contains all routines and parameters for the application of
!!              averaging techniques and the providing of mean elements.
!!
!!              \par Includes
!!
!!

!!
 neptune_error_handling.f90

!!
 slam_math.f90

!!
 slam_astro_conversions.f90

!!
 slam_orbit_types.f90

!!
 slam_types.f90

!!
 slam_units.f90

!!

!!
!> @copyright   Institute of Space Systems / TU Braunschweig
!!
!> @anchor      averaging
!!
!!------------------------------------------------------------------------------------------------
module averaging

  use slam_astro_conversions,  only: getOrbitType
  use slam_error_handling,     only: FATAL, E_KEPLER_UNITS, isControlled, hasToReturn, checkIn, checkOut
  use neptune_error_handling,  only: setNeptuneError, E_MEAN_ELS_MISSING
  use slam_math,               only: deg2rad
  use slam_orbit_types,        only: kepler_t, assignment(=)
  use slam_types,              only: dp
  use slam_units,              only: UNIT_DEG, UNIT_RAD, UNIT_KM, UNIT_M

  implicit none

  private

  logical :: isAvailable = .false.    ! flag for the availability of a current mean state (kepler elements)

  type(kepler_t) :: meanState

! Averaging is not a class yet, but meanState and isAvailable should be private per thread
!$omp threadprivate(isAvailable,meanState)

  !===================================
  !
  ! public methods
  !
  !------------------
  public :: getMeanElements
  public :: setMeanElements

contains

!>-----------------------------------------------------------------------------------------------
!> @brief       Flag for the availability of mean elements as global variable of this module
!!
!> @author      Vitali Braun
!!
!> @return      .true. when mean elements are available
!!
!> @date

!!
 30.07.2013 (initial design)

!!

!!
!> @anchor      isAvailableMeanState
!!
!!------------------------------------------------------------------------------------------------
  logical function isAvailableMeanState()

    isAvailableMeanState = isAvailable

  end function isAvailableMeanState

!>-----------------------------------------------------------------------------------------------
!> @brief       Returns the current mean elements
!!
!> @author      Vitali Braun
!!
!> @return      mean keplerian elements as kepler_t
!!
!! @date

!!
 30.07.2013 (initial design)

!!

!!
!! @anchor      getMeanElements
!!
!!------------------------------------------------------------------------------------------------
  type(kepler_t) function getMeanElements()

    character(len=*), parameter :: csubid = "getMeanElements"
    real(dp) :: zero = 0.0_dp

    getMeanElements = zero

    if(isControlled()) then
      if(hasToReturn()) return
      call checkIn(csubid)
    end if

    if(isAvailable) then
      getMeanElements = meanState

    else
      call setNeptuneError(E_MEAN_ELS_MISSING, FATAL)
      getMeanElements = zero
      return

    end if

    !** done!
    if(isControlled()) then
      call checkOut(csubid)
    end if

    return

  end function getMeanElements

!>-----------------------------------------------------------------------------------------------
!> @brief       Set mean elements
!!
!> @author      Vitali Braun
!!
!> @date

!!
 30.07.2013 (initial design)

!!

!!
!> @param[in]   kep     Kepler elements
!!
!> @anchor      setMeanElements
!!
!!------------------------------------------------------------------------------------------------
  subroutine setMeanElements(kep)

    type(kepler_t), intent(in) :: kep

    character(len=*), parameter :: csubid = 'setMeanElements'
    integer :: orbtype

    if(isControlled()) then
      if(hasToReturn()) return
      call checkIn(csubid)
    end if

    meanState = kep

    !** check for RAD/DEG units
    select case(meanState%angles_unit)

      case(UNIT_DEG) ! convert to radian
        meanState%inc     = meanState%inc*deg2rad       ! inclination
        meanState%raan    = meanState%raan*deg2rad      ! right ascension of ascending node
        meanState%aop     = meanState%aop*deg2rad       ! argument of pericenter
        meanState%man     = meanState%man*deg2rad       ! mean anomaly
        meanState%ecan    = meanState%ecan*deg2rad      ! eccentric anomaly
        meanState%tran    = meanState%tran*deg2rad      ! true anomaly
        meanState%arglat  = meanState%arglat*deg2rad    ! argument of true latitude
        meanState%lonper  = meanState%lonper*deg2rad    ! longitude of perigee
        meanState%truelon = meanState%truelon*deg2rad   ! true longitude

        meanState%angles_unit = UNIT_RAD

      case(UNIT_RAD)
        ! nothing happens

      case default
        call setNeptuneError(E_KEPLER_UNITS, FATAL)
        return

    end select

    !** check for SMA unit
    select case(meanState%sma_unit)

      case(UNIT_KM)
        ! nothing happens..

      case(UNIT_M)
        meanState%sma = meanState%sma/1.d3
        meanState%sma_unit = UNIT_KM

      case default
        call setNeptuneError(E_KEPLER_UNITS, FATAL)
        return

    end select

    !** determine orbit type and compute missing values
    call getOrbitType(meanState,orbtype)

    !** done!
    if(isControlled()) then
      call checkOut(csubid)
    end if

    isAvailable = .true.
    return

  end subroutine setMeanElements

end module averaging

max_stars_repo_name: mmoeckel/neptune

max_issues_repo_name: mmoeckel/neptune

max_forks_repo_name: mmoeckel/neptune


text:
C SUBROUTINE CMPINV           (NEXT)
      SUBROUTINE CMPINV (A, B, C, D, N)
C
C     CMPINV - COMPLEX MATRIX OPERATION THAT FINDS INVERSE, C+ID, OF
C        A+IB WHERE A,B,C, AND D ARE REAL SQUARE MATRICES.  CMPINV
C        CALLS MATINV TO INVERT REAL MATRICES AND SQMULT TO MULTIPLY
C        SQUARE MATRICES IN ORDER TO COMPUTE THE COMPLEX INVERSE
C                      C=INV(A + B*(INV(A))*B)
C                      D=-C*B*(INV(A))
C        A, B, C, AND D ARE ALL REAL N BY N MATRICIES.
C        IF A IS SINGULAR, EXCHANGE A AND B (FACTOR OUT I) AND TRY AGAIN
C     THIS ROUTINE IS SHARE NO. ANL F103 BY B.S.GARBOW OF ARGONNE
C     NATIONAL LABORATORY, ARGONNE, ILLINOIS
C
      DIMENSION A (20, 20), B (20, 20), C (20, 20), D (20, 20)
C
      COMMON /PON /INDEX (20, 2), IPIVOT (20), PIVOT (20)
C
C.....INVERT -A.........................................................
      M = N
      L = 1
 100  IF (M .LT. 1) GO TO 160
      DO 105 I = 1, M
      IF (M .LT. 1) GO TO 165
      DO 105 J = 1, M
 105  D (I, J) = - A (I, J)
 165  CONTINUE
 160  CONTINUE
      CALL MATINV (D, M)
C.....CHECK IF INV(-A) EXISTS (IF NOT GO TO 6)..........................
      IF (M .LT. 1) GO TO 170
      DO 110 K = 1, M
      IF (IPIVOT (K) .NE. 1) GO TO 125
 110  CONTINUE
 170  CONTINUE
C.....COMPUTE C=INV(A - B*(INV(-A))*B)..................................
      CALL SQMULT (B, D, D, M)
      CALL SQMULT (D, B, C, M)
      IF (M .LT. 1) GO TO 175
      DO 115 I1 = 1, M
      IF (M .LT. 1) GO TO 180
      DO 115 J1 = 1, M
 115  C (I1, J1) = A (I1, J1) - C (I1, J1)
 180  CONTINUE
 175  CONTINUE
      CALL MATINV (C, M)
C.....CHECK IF C EXISTS (IF NOT GO TO 300)..............................
      IF (M .LT. 1) GO TO 185
      DO 120 K1 = 1, M
      IF (IPIVOT (K1) .NE. 1) GO TO 150
 120  CONTINUE
 185  CONTINUE
C.....COMPUTE D=C*B*(INV(-A))...........................................
      CALL SQMULT (C, D, D, M)
      IF(L.NE.2) GO TO 145
      GO TO 135
C.....A IS SINGULAR, INTERCHANGE A AND B AND TRY AGAIN..................
 125  IF (M .LT. 1) GO TO 190
      DO 130 I2 = 1, M
      IF (M .LT. 1) GO TO 195
      DO 130 J2 = 1, M
      TEMP = A (I2, J2)
      A (I2, J2) = B (I2, J2)
 130  B (I2, J2) = TEMP
 195  CONTINUE
 190  CONTINUE
C.....IF A AND B BOTH SINGULAR, GO TO 310...............................
      IF (L .EQ. 2) GO TO 155
      L = 2
      GO TO 100
C.....RE-INTERCHANGE A AND B, AND INTERCHANGE C AND D WITH SIGNS CHANGED
C     DUE TO THE -I FACTORED OUT WHEN A AND B ORIGINALLY INTERCHANGED...
 135  IF (M .LT. 1) GO TO 200
      DO 140 I3 = 1, M
      IF (M .LT. 1) GO TO 205
      DO 140 J3 = 1, M
      TEMP = A (I3, J3)
      A (I3, J3) = B (I3, J3)
      B (I3, J3) = TEMP
      TEMP = - C (I3, J3)
      C (I3, J3) = - D (I3, J3)
 140  D (I3, J3) = TEMP
 205  CONTINUE
 200  CONTINUE
 145  RETURN
C
C.....ERROR STOP - NO INVERSE EXISTS....................................
 150  WRITE (* , 1500)
1500  FORMAT ('0MATRIX C IN SUBROUTINE CMPINV DOES NOT EXIST.'//)
      STOP
C.....ERROR STOP - INVERSE CANNOT BE FOUND BY THIS METHOD...............
 155  WRITE (* , 1502)
1502  FORMAT ('0MATRICES A AND B BOTH SINGULAR IN SUBR CMPINV.'//)
      STOP
      END

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
test1/test.cpp
test1/tb.cpp
test1/main.cpp


! { dg-do compile }
! Tests the fix for PR31550 in which pointers to derived type components
! were being TREE-SSA declared in the wrong order and so in the incorrect
! context.
!
! Contributed by Daniel Franke <dfranke@gcc.gnu.org>
!
MODULE class_dummy_atom_types
TYPE :: dummy_atom_list
  TYPE(dummy_atom), DIMENSION(:), POINTER :: table
  INTEGER                                 :: nused
END TYPE

TYPE :: dummy_atom
  TYPE(dummy_atom_private), POINTER :: p
END TYPE

TYPE :: dummy_atom_private
  TYPE(dummy_atom_list)       :: neighbors
END TYPE
END MODULE

MODULE class_dummy_atom_list
USE class_dummy_atom_types, ONLY: dummy_atom_list

INTERFACE
  SUBROUTINE dummy_atom_list_init_copy(this, other)
    USE class_dummy_atom_types, ONLY: dummy_atom_list
    TYPE(dummy_atom_list), INTENT(out) :: this
    TYPE(dummy_atom_list), INTENT(in)  :: other
  END SUBROUTINE
END INTERFACE

INTERFACE
  SUBROUTINE dummy_atom_list_merge(this, other)
    USE class_dummy_atom_types, ONLY: dummy_atom_list
    TYPE(dummy_atom_list), INTENT(inout) :: this
    TYPE(dummy_atom_list), INTENT(in)    :: other
  END SUBROUTINE
END INTERFACE
END MODULE

SUBROUTINE dummy_atom_list_init_copy(this, other)
  USE class_dummy_atom_list, ONLY: dummy_atom_list, dummy_atom_list_merge

  TYPE(dummy_atom_list), INTENT(out) :: this
  TYPE(dummy_atom_list), INTENT(in)  :: other

  this%table(1:this%nused) = other%table(1:other%nused)
END SUBROUTINE


module tools
  use types
  ! ----------------------------------------------------------------
  ! tools:
  !   This module contains all the functions and subroutines
  ! author:
  !   Exam-roll -- 20419PHY083
  ! env details:
  !   compiler -- gfortran v10.2.0
  !   OS       -- macOS 12.0 x86_64
  !   Kernel   -- Darwin 21.0.0
  ! Note:
  !   - depends on the module "types"
  ! ----------------------------------------------------------------
  implicit none

contains
  function fn(x) result(z)
    ! ----------------------------------------------------------------
    ! Single variable function:
    !   function with one variable that needs to be calculated
    ! input:
    !   x -- input of the function
    ! output:
    !   z -- output of the function f(x)
    ! ----------------------------------------------------------------
    implicit none
    real(dp), intent(in) :: x
    real(dp) :: z

    z = 2.0*x**3.0 - 3.0*x**2.0 + 4.0*x - 5.0
    ! z = 0.5*exp(-((1 + x)/2.0)**2.0)
    ! z = 1/(x + 2)
  end function fn

  function f(x, y) result(z)
    ! ----------------------------------------------------------------
    ! Two variable function:
    !   function with two variables that needs to be calculated
    ! inputs:
    !   x,y -- inputs of the function
    ! output:
    !   z   -- output of the function f(x,y)
    ! ----------------------------------------------------------------
    implicit none
    real(dp), intent(in) :: x, y
    real(dp) :: z

    ! z = exp(x) + y
    z = -x*y
    ! z = sqrt(1 - x**2.0)
  end function f

  subroutine oela(x0, y0, x_point, h)
    ! ----------------------------------------------------------------
    ! Euler method:
    !   Finds the value of y at a given x by using a picewise linear
    !   approximation of the solution
    ! inputs:
    !   h         -- interval size (keep it small for better result)
    !   x0        -- initial value of x (also used as output)
    !   y0        -- initial value of y (also used as output)
    !   x_point   -- the value of x at which y (f(x)) is to be calculated
    ! internals:
    !   x_new     -- new value of x's calculated
    !   y_new     -- new value of y's calculated
    !   k1        -- slope
    ! outputs:
    !   y0, x0    -- values of the required y at given x
    ! ----------------------------------------------------------------
    implicit none
    real(dp), intent(in) :: x_point, h
    real(dp) :: x0, y0
    real(dp) :: k1, x_new, y_new
    integer :: itn

    x_new = x0
    y_new = y0
    itn = 0

    do while (x_new <= x_point)
      ! displaying the values
      write (*, form_table3) itn, x_new, y_new
      x0 = x_new
      y0 = y_new

      ! valus of slope
      k1 = f(x_new, y_new)
      ! calculating the y
      y_new = y_new + h*k1

      ! preparing for the next value
      x_new = x_new + h
      itn = itn + 1
    end do
  end subroutine oela

  subroutine mod_oela(x0, y0, x_point, h)
    ! ----------------------------------------------------------------
    ! Modified Euler method:
    !   Finds the value of y at a given x also known as Heun's method
    ! inputs:
    !   h         -- interval size (keep it small for better result)
    !   x0        -- initial value of x (also used as output)
    !   y0        -- initial value of y (also used as output)
    !   x_point   -- the value of x at which y (f(x)) is to be calculated
    ! internals:
    !   x_new     -- new value of x's calculated
    !   y_new     -- new value of y's calculated
    !   k1, k2    -- slope
    ! outputs:
    !   y0, x0    -- values of the required y at given x
    ! ----------------------------------------------------------------
    implicit none
    real(dp), intent(in) :: h, x_point
    real(dp) :: x0, y0
    real(dp) :: k1, k2, k, x_new, y_new
    integer :: itn

    x_new = x0
    y_new = y0
    itn = 0

    do while (x_new <= x_point)
      ! displaying the values
      write (*, form_table3) itn, x_new, y_new
      x0 = x_new
      y0 = y_new

      ! values of slope
      k1 = f(x_new, y_new)
      k2 = f(x_new + h, y_new + h*k1)

      ! weighted average of the slopes
      k = (k1 + k2)/2.0

      ! calculating the y
      y_new = y_new + h*k

      ! preparing for the next value
      x_new = x_new + h
      itn = itn + 1
    end do
  end subroutine mod_oela

  subroutine g_quad(k, a, b, int_val)
    ! ----------------------------------------------------------------
    ! Gaussian quadrature method:
    !   Integrates a given function by fitting a polynomial at unequal
    !   intervals whose abscissa and weights are calculated using the
    !   zeros of the legendre polynomial and Newton's method.
    ! inputs:
    !   k         -- number of points needed
    !   x()       -- the arrays of x's
    !   y()       -- the arrays of corresponding y's (f(x)'s)
    !   x_point   -- the value of x at which f'(x) is to be calculated
    ! internal:
    !   h         -- class height
    !   loc       -- determines the position of x_point on the table
    !                and gets the respective index for calculation
    !   y         -- value of the integral before scaling
    ! output:
    !   int_val   -- the result of the derivation after scaling
    ! Note:
    !   - used function fn() because only 1 variable
    !   - used subroutine absc_wght() for abscissa and weights
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: k
    integer, intent(in) :: a, b
    real(dp), intent(out) :: int_val

    real(dp), dimension(k) :: x, w
    real(dp) :: fun, y, p, q
    integer :: i

    call absc_wght(k, x, w) ! calculating the weighting factors

    ! printing the above weighting factors
    write (*, form_header2) 'abscissa', 'weights'
    write (*, line2) '--------', '--------'
    do i = 1, k
      write (*, form_table2) x(i), w(i)
    end do

    ! main integration
    y = 0.0
    p = (a + b)/2.0; q = (b - a)/2.0

    ! loop to calculate the values for each weighted factor
    do i = 1, k
      fun = fn(p + q*x(i))
      y = y + w(i)*fun
    end do

    int_val = q*y ! final result
  end subroutine g_quad

  subroutine absc_wght(points, abscissa, weights)
    ! ----------------------------------------------------------------
    ! Gaussian quadrature Weighting Factors:
    !   calculate the abscissa and weights  for the gaussian quadrature
    !   method {g_quad()} using the zeros of the legendre polynomial
    ! inputs:
    !   points    -- number of data points
    !   abscissa  -- the abscissa
    !   weights   -- the weights
    ! internals:
    !   p0        -- Pn coefficient of the legendre polynomial
    !   p1        -- Pn-1 coefficient of the legendre polynomial
    !   p2        -- Pn-2 coefficient of the legendre polynomial
    !   m         -- finds the mid-point
    !   x         -- initial guess of the abscissa which is further refined
    !               using the Newton-rhapson method
    !   p0_dash   -- the derivative of the legendre polynomial
    ! Note:
    !   - Newton's method is used to refine the value of abscissa (x)
    !   - since the weighted factors are symmetric we can break it into half
    !     and calculate only one half and then set the other half
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: points
    real(dp), dimension(points), intent(out) :: abscissa, weights

    real(dp), parameter :: eps = 1.0e-15
    integer :: i, j, m
    real(dp) :: p0, p1, p2, p0_dash, x, x_old

    if (mod(points, 2) == 0) then ! finding the mid-point
      m = points/2
    else
      m = (points + 1)/2
    end if

    ! loop to calculate the zeros of legendre polynomial
    do i = 1, m
      x = cos(pi*(i - 0.25)/(points + 0.5)) ! initial guess of x

      do while (abs(x - x_old) > eps)
        p0 = 1.0
        p1 = x

        ! loop to calculate the legendre polynomial
        do j = 1, points
          p2 = p1
          p1 = p0
          p0 = ((2.0*j - 1.0)*x*p1 - (j - 1.0)*p2)/j ! legendre polynomial
        end do

        p0_dash = points*(x*p0 - p1)/(x**2.0 - 1.0) ! legendre polynomial derivative

        x_old = x
        x = x_old - p0/p0_dash ! newton's method for refining x

        ! arranging the abscissa
        abscissa(i) = -x ! the -abscissa's
        abscissa(points + 1 - i) = x ! the +abscissa's

        ! arranging the weights
        weights(i) = (2.0*(1 - x**2.0))/(points*p1)**2.0 ! the weights for -abscissa
        weights(points + 1 - i) = weights(i) ! the weights for +abscissa
      end do
    end do
  end subroutine absc_wght

  subroutine rungkut_2(x0, y0, x_point, h)
    ! ----------------------------------------------------------------
    ! Runge-Kutte two point method:
    !   Finds the value of y at a given x also known as Polygon method
    ! inputs:
    !   h         -- interval size (keep it small for better result)
    !   x0        -- initial value of x (also used as output)
    !   y0        -- initial value of y (also used as output)
    !   x_point   -- the value of x at which y (f(x)) is to be calculated
    ! internals:
    !   x_new     -- new value of x's calculated
    !   y_new     -- new value of y's calculated
    !   k1, k2    -- slopes
    ! outputs:
    !   y0, x0    -- values of the required y at given x
    ! ----------------------------------------------------------------
    implicit none
    real(dp), intent(in) :: h, x_point
    real(dp) :: x0, y0
    real(dp) :: k1, k2, x_new, y_new
    integer :: itn

    x_new = x0
    y_new = y0
    itn = 0

    do while (x_new <= x_point)
      ! displaying the values
      write (*, form_table3) itn, x_new, y_new
      x0 = x_new
      y0 = y_new
      ! value of slopes
      k1 = f(x_new, y_new)
      k2 = f(x_new + h/2.0, y_new + k1*h/2.0)
      ! calculating the y
      y_new = y_new + h*k2
      ! preparing for next values
      x_new = x_new + h
      itn = itn + 1
    end do
  end subroutine rungkut_2

  subroutine rungkut_4(x0, y0, x_point, h)
    ! ----------------------------------------------------------------
    ! Runge-Kutte two point method:
    !   Finds the value of y at a given x also known as Polygon method
    ! inputs:
    !   h         -- interval size (keep it small for better result)
    !   x0        -- initial value of x (also used as output)
    !   y0        -- initial value of y (also used as output)
    !   x_point   -- the value of x at which y (f(x)) is to be calculated
    ! internals:
    !   x_new     -- new value of x's calculated
    !   y_new     -- new value of y's calculated
    !   k1,k2,k3,k4  -- slopes
    ! outputs:
    !   y0, x0    -- values of the required y at given x
    ! ----------------------------------------------------------------
    implicit none
    real(dp), intent(in) :: x_point, h
    real(dp) :: x0, y0
    real(dp) :: k1, k2, k3, k4, k, x_new, y_new
    integer :: itn

    x_new = x0
    y_new = y0
    itn = 0

    do while (x_new <= x_point)
      ! displaying the values
      write (*, form_table3) itn, x_new, y_new
      x0 = x_new
      y0 = y_new

      ! values of slopes
      k1 = f(x_new, y_new)
      k2 = f(x_new + h/2.0, y_new + k1*h/2.0)
      k3 = f(x_new + h/2.0, y_new + k2*h/2.0)
      k4 = f(x_new + h, y_new + k3*h)

      ! weighted average of the slopes
      k = (k1 + 2.0*(k2 + k3) + k4)/6.0

      ! calculating the y
      y_new = y_new + k*h

      ! preparingg for the next value
      x_new = x_new + h
      itn = itn + 1
    end do
  end subroutine rungkut_4

  subroutine neum_diff(n, x, y, x_point, poly, deriv)
    ! ----------------------------------------------------------------
    ! Newton's differentiation:
    !   Calculates the differentiation of y at a given x using newton's
    !   interpolation polynomial and the difference table calculated
    !   by the subroutine newinter()
    ! inputs:
    !   n         -- number of points for interpolation
    !   x()       -- the arrays of x's
    !   y()       -- the arrays of corresponding y's (f(x)'s)
    !   x_point   -- the value of x at which f'(x) is to be calculated
    ! internal:
    !   m         -- the mid-point of the table of n points
    !   h         -- class height
    !   loc       -- determines the position of x_point on the table
    !               and gets the respective index for calculation
    ! output:
    !   deriv     -- the result of the derivation
    !   poly      -- the used interpolation polynomial
    ! Note:
    !   - used the subroutine newinter
    !   - used forward difference table and modified the formula a bit
    !     to calculate the backward differentiation.
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: n
    real(dp), dimension(n), intent(in) :: x, y
    real(dp), intent(in) :: x_point
    real(dp), intent(out) :: deriv
    character(len=8), intent(out) :: poly

    real(dp), dimension(n, n) :: del_f
    real(dp) :: h, sum_1, sum_2, x_sum
    integer :: a, b, m, loc

    ! condition to get the mid-point of the data points
    if (mod(n, 2) == 0) then
      m = n/2
    else
      m = (n + 1)/2
    end if

    loc = minloc(x, mask=(x >= x_point), dim=1) ! position of an x in the array

    ! check if x_point is outside the given data points
    if ((loc == 0) .or. (loc == n + 1)) then
      stop 'your {x_point} is outside the given data points'
    end if

    call newinter(n, x, y, del_f) ! calculating the difference table

    ! conditon for forward and backward differentiation
    if (x_point <= x(m)) then ! forward case
      poly = 'forward'
      if (loc == 1) then
        a = loc
        b = a + 1
      else
        a = loc
        b = a - 1
      end if

      ! calculation of the derivative
      h = abs(x(a) - x(b))
      x_sum = ((x_point - x(b)) + (x_point - x(a)))
      sum_1 = del_f(b, 1)/h
      sum_2 = del_f(b, 2)/(2.0*h*h)

    else if (x_point > x(m)) then ! backward case
      poly = 'backward'
      if (loc == n) then
        a = loc
        b = a - 1
      else
        a = loc
        b = a + 1
      end if

      ! calculation of the derivative
      h = abs(x(a) - x(b))
      x_sum = ((x_point - x(b)) + (x_point - x(a)))
      sum_1 = del_f(b, 1)/h
      sum_2 = del_f(b - 1, 2)/(2.0*h*h)
    end if

    deriv = sum_1 + sum_2*x_sum ! final result
  end subroutine neum_diff

  subroutine newinter(n, x, y, del_f)
    ! ----------------------------------------------------------------
    ! Newton's interpolation table:
    !   Calculates the forward differences and gives us the F.D. table
    ! inputs:
    !   n         -- number of points for interpolation
    !   x()       -- the arrays of x's (used just for tabulation)
    !   y()       -- the arrays of corresponding y's (f(x)'s)
    ! output:
    !   del_f(,)  -- the forward difference coefficients
    ! Note:
    !   - This same table is used for backward difference because it is
    !     same as forward difference, just the position is different.
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: n
    real(dp), dimension(n), intent(in) :: y, x
    real(dp), dimension(n, n), intent(out) :: del_f

    real(dp), dimension(n, n) :: del_y
    integer :: i, j, p

    ! loop to calculate the differences
    do j = 1, n
      do i = 1, n - j
        if (j == 1) then ! calculating the 1st difference column
          del_y(i, 1) = y(i + 1) - y(i)
        else ! calculate the other coloumns
          del_y(i, j) = del_y(i + 1, j - 1) - del_y(i, j - 1)
        end if
      end do
    end do

    del_f = del_y

    ! printing the difference table
    print *, 'the difference table:'
    write (*, form_header3n) 'x', 'y', "del_f's"
    write (*, line3n) '------', '------', '-------------------------'
    do i = 1, n
      write (*, form_tableip) x(i), y(i), (del_f(i, p), p=1, n - i)
    end do
  end subroutine newinter

  subroutine jac(n, a, x, b, err, itn)
    ! ----------------------------------------------------------------
    ! Jacobi method:
    !   calculate the single value of x,y and z using the iterative
    !   method of simultaneous displacememnts
    ! inputs:
    !   n         -- number of data points (dimension of the matrix)
    !   a         -- matrix of the coefficients A
    !   b         -- matrix of the constants B
    !   x(:,1)    -- initial x1,x2,x3
    !   itn       -- iteration number
    ! internals:
    !   sum_ax()  -- sum of the off diagonal terms of A
    !   aberror   -- difference of the values of old x's with new x's
    ! output:
    !   err       -- the eucledian norm of "error"
    !   x(:,2)    -- final x1,x2,x3
    ! Note:
    !   - use the subroutine dominance() to check if the solution will
    !     converge or not
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: n, itn
    real(dp), dimension(n), intent(in) :: b
    real(dp), dimension(n, n), intent(in) :: a
    real(dp), intent(out) :: err
    real(dp), dimension(n, 2) :: x

    real(dp), dimension(n) :: sum_ax, aberror
    integer :: p, q

    !loop to calculate x1,x2,x3
    do p = 1, n
      sum_ax(p) = 0.0
      do q = 1, n
        if (q /= p) then
          sum_ax(p) = sum_ax(p) + a(p, q)*x(q, 1)
        end if
      end do
      x(p, 2) = (b(p) - sum_ax(p))/a(p, p)
    end do

    ! error calculation
    aberror = x(:, 1) - x(:, 2)
    err = norm2(aberror) ! taking the eucledian norm

    ! printing the values
    write (*, form_table6) itn, (x(p, 2), p=1, n), err
  end subroutine jac

  subroutine dominance(n, a)
    ! ----------------------------------------------------------------
    ! Diagonally dominance checker:
    !   checks if the matrix is strictly diagonally dominant
    ! inputs:
    !   n         -- dimension of the matrix
    !   a         -- matrix of the coefficients A
    ! internals:
    !   row_sum   -- sum of the off diagonal row elements
    ! Note:
    !   - used for the jacobi method equation solver.
    !   - if dominance test failed, need to make the matrix dominant
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: n
    real(dp), dimension(n, n), intent(in) :: a

    real(dp) :: row_sum
    integer :: p, q

    ! loop to calculate the sum of off diagonal terms
    do p = 1, n
      row_sum = 0.0
      do q = 1, n
        if (q /= p) then
          row_sum = row_sum + abs(a(p, q)) ! off-diag element sum
        end if
      end do

      ! condition for strictly diagonally dominant
      if (abs(a(p, p)) < row_sum) then
        stop "-- !! the matirx is NOT strictly diagonally dominant !! --"
      end if
    end do
  end subroutine dominance

  subroutine monte_pi(sq_points, points_inside, approx_pi)
    ! ----------------------------------------------------------------
    ! Monte Carlo Pi:
    !   Calculates the approximate value of pi using Monte Carlo method
    ! inputs:
    !   sq_points     -- total points inside square including the circle
    !   points_inside -- total points inside the circle
    ! internals:
    !   rand_x, rand_y -- uniform random number in the range [0,1)
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: sq_points
    integer, intent(out) :: points_inside
    real(dp), intent(out) :: approx_pi

    integer :: i, cr_points
    real(dp) :: x, y, rand_x, rand_y, cr_radi

    call random_seed() ! to maintain the state random_number generator

    do i = 1, sq_points
      ! generating uniform random numbers in the range [0,1)
      call random_number(rand_x)
      call random_number(rand_y)

      ! scaling in the range [a,b]: x = a + (b-a)*x
      x = -1 + 2*rand_x
      y = -1 + 2*rand_y

      cr_radi = x**2.0 + y**2.0 ! radius of the circle

      ! condition to check if the hits were perfect
      if (cr_radi <= 1.0) then
        cr_points = cr_points + 1
      end if
    end do

    ! final result
    points_inside = cr_points
    approx_pi = 4.0*(float(cr_points)/float(sq_points))
  end subroutine monte_pi

  subroutine monte_int(a, b, s_size, approx_int)
    ! ----------------------------------------------------------------
    ! Monte Carlo Integration:
    !   Calculates the approximate value of integral using MC method
    ! inputs:
    !   a,b           -- lower and upper limit
    !   s_size        -- size of the sample
    ! internals:
    !   rand_x, rand_y -- uniform random number in the range [0,1)
    !   sum_f         -- sum of the function from 1 to sample size
    ! output:
    !   approx_int    -- the approximate value of the integral
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: a, b, s_size
    real(dp), intent(out) :: approx_int

    real(dp) :: x, y, sum_f, rand_x, rand_y
    integer :: i

    call random_seed() ! to maintain the state random_number generator

    ! loop to calculate the integral
    do i = 1, s_size
      ! generating uniform random numbers in the range [0,1)
      call random_number(rand_x)
      call random_number(rand_y)

      ! scaling in the range [a,b]: x = a + (b-a)*x
      x = a + rand_x*(b - a)
      y = a + rand_y*(b - a)

      sum_f = sum_f + f(x, y) ! summing all the values of the given function
    end do

    approx_int = float(b - a)*(sum_f/float(s_size)) ! final result of integral
  end subroutine monte_int

  subroutine garphmonte_pi(sq_points, points_inside, approx_pi)
    ! ----------------------------------------------------------------
    ! Monte Carlo Pi:
    !   Calculates the approximate value of pi using Monte Carlo method
    ! inputs:
    !   sq_points     -- total points inside square including the circle
    !   points_inside -- total points inside the circle
    ! internals:
    !   rand_x, rand_y -- uniform random number in the range [0,1)
    ! Note:
    !   - this subroutine can be used to get the points in a output file
    !     and can be plotted using a graphing tool
    ! ----------------------------------------------------------------
    implicit none
    integer, intent(in) :: sq_points
    integer, intent(out) :: points_inside
    real(dp), intent(out) :: approx_pi

    integer :: i, cr_points
    real(dp) :: cr_radi
    real(dp) :: x, y, rand_x, rand_y

    call random_seed() ! to maintain the state random_number generator

    open (1, file='values_inside.dat')
    open (2, file='values_outside.dat')

    do i = 1, sq_points
      ! generating uniform random numbers in the range [0,1)
      call random_number(rand_x)
      call random_number(rand_y)

      ! scaling in the range [a,b]: x = a + (b-a)*x
      x = -1 + 2*rand_x
      y = -1 + 2*rand_y

      cr_radi = x**2.0 + y**2.0 ! radius of the circle

      ! condition to check if the hits were perfect
      if (cr_radi <= 1.0) then
        cr_points = cr_points + 1
        write (1, "(2f15.5)") x, y
      else
        write (2, "(2f15.5)") x, y
      end if
    end do

    close (1); close (2)

    ! final result
    points_inside = cr_points
    approx_pi = 4.0*(float(cr_points)/float(sq_points))
  end subroutine garphmonte_pi

end module tools