text:
C$Procedure  ZZEKJSQZ ( Private: EK, join row set squeeze )

      SUBROUTINE ZZEKJSQZ ( JRSBAS )

C$ Abstract
C
C     SPICE Private routine intended solely for the support of SPICE
C     routines.  Users should not call this routine directly due
C     to the volatile nature of this routine.
C
C     Compress a join row set by eliminating segment vectors for
C     which there are no corresponding row vectors.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     EK
C
C$ Keywords
C
C     EK
C     PRIVATE
C
C$ Declarations

      INCLUDE 'ekjrs.inc'
      INCLUDE 'ekqlimit.inc'

      INTEGER               JRSBAS

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     JRSBAS     I   Scratch area base address of join row set.
C
C$ Detailed_Input
C
C     JRSBAS         is the base address, in the scratch area, of a
C                    join row set to be compressed.
C
C$ Detailed_Output
C
C     None.  See $Particulars for a description of the effect of this
C     routine.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  If JRSBAS is not the base address of a structurally valid
C         join row set, the results of this routine will be
C         unpredictable.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine operates by side effects:  it modifies the join row
C     set designated by the input argument JRSBAS.  Every row vector
C     marked for deletion is removed.  Every empty segment vector is
C     removed, along with the row count and row vector base for that
C     segment vector.  The join row set is compressed to remove all
C     gaps.   All counts are updated to reflect the updated join row
C     set.
C
C     The purpose of the compression performed by this routine is to
C     save work during joins by reducing the size of the cartesian
C     products of sets of segment vectors.  Also, special cases
C     involving null segment vectors can be avoided by this clean-up
C     mechanism.  Finally, it may be possible to save space in the EK
C     scratch area freed by the compression.
C
C$ Examples
C
C     See EKSRCH.
C
C$ Restrictions
C
C     1) Relies on the EK scratch area.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman       (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.1.0, 07-AUG-2006 (NJB)
C
C        Bug fix:  added intialization of variable NRVDEL to support
C                  operation under the Macintosh Intel Fortran
C                  compiler. Note that this bug did not affect
C                  operation of this routine on other platforms.
C
C-    SPICELIB Version 1.0.0, 10-OCT-1995 (NJB)
C
C-&

C$ Revisions
C
C-    SPICELIB Version 1.1.0, 07-AUG-2006 (NJB)
C
C        Bug fix:  added intialization of variable NRVDEL to support
C                  operation under the Macintosh Intel Fortran
C                  compiler. Note that this bug did not affect
C                  operation of this routine on other platforms. The
C                  statement referencing the uninitialized variable
C                  was:
C
C           IF (  ( RC .EQ. 0 ) .OR. ( NRVDEL .EQ. RC )  ) THEN
C
C        In the previous version of the code, NRVDEL is uninitialized
C        when NRVDEL is 0.  NRVDEL *is* initialized when RC is
C        non-zero, so the logical value of the IF expression is not
C        affected by the lack of proper intialization.
C
C        However, the Intel Fortran compiler for the Mac flags a runtime
C        error when the above code is exercised.  So NRVDEL is now
C        initialized prior to the above IF statement.
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               FAILED

C
C     Local variables
C
      INTEGER               CNTLOC
      INTEGER               DELTA
      INTEGER               I
      INTEGER               J
      INTEGER               NSVDEL
      INTEGER               NRVDEL
      INTEGER               NEWNSV
      INTEGER               NR
      INTEGER               NRLOC
      INTEGER               NSV
      INTEGER               NSVLOC
      INTEGER               NTAB
      INTEGER               NTLOC
      INTEGER               PCPAIR ( 2 )
      INTEGER               PTARG
      INTEGER               PTBASE
      INTEGER               PTRLOC
      INTEGER               RC
      INTEGER               RBASE
      INTEGER               ROWVEC ( MAXTAB+1 )
      INTEGER               RTARG
      INTEGER               RVSIZE
      INTEGER               SEGVEC ( MAXTAB )
      INTEGER               SETBAS
      INTEGER               SIZE
      INTEGER               SIZLOC
      INTEGER               SVBASE
      INTEGER               SVSIZE
      INTEGER               VTARG

C
C     Use discovery check-in.
C
C
C     Look up the counts that are of interest:
C
C       -- The table count
C       -- The segment vector count
C       -- The join row set size
C
C     Save the address of each count.
C
      SIZLOC  =  JRSBAS + JSZIDX
      NSVLOC  =  JRSBAS + JSCIDX
      NTLOC   =  JRSBAS + JTCIDX

      CALL ZZEKSRD ( SIZLOC, SIZLOC, SIZE   )
      CALL ZZEKSRD ( NTLOC,  NTLOC,  NTAB   )
      CALL ZZEKSRD ( NSVLOC, NSVLOC, NSV    )

      IF ( FAILED() ) THEN
         RETURN
      END IF

C
C     Set the sizes of segment and row vectors.
C
      SVSIZE  =  NTAB
      RVSIZE  =  NTAB + 1

C
C     For each segment vector, obtain the row count.  Clean up after
C     null segment vectors:  compress out the space allocated for their
C     row vector pointers.  Keep track of the number of deletions.
C
      NSVDEL  =  0
      NRVDEL  =  0
      VTARG   =  JRSBAS + JSVBAS

      DO I = 1, NSV
C
C        The location of the row count is CNTLOC.  The row vector base
C        pointer precedes the row count.
C
         CNTLOC  =  JRSBAS + JSVBAS + NSV*SVSIZE + 2*(I-1) + 2
         PTRLOC  =  CNTLOC - 1

         CALL ZZEKSRD ( CNTLOC, CNTLOC, RC )

         IF ( RC .GT. 0 ) THEN
C
C           The row vector set for this segment vector is non-empty.
C           scan the rows, looking for those marked for deletion, and
C           update the row count to reflect the number of rows that
C           we're going to keep.
C
            CALL ZZEKSRD ( PTRLOC, PTRLOC, SETBAS )

            NRVDEL = 0

            DO J = 1, RC

               RBASE  =  JRSBAS + SETBAS + (J-1)*RVSIZE

               CALL ZZEKSRD ( RBASE+1, RBASE+1, ROWVEC(1) )

               IF ( ROWVEC(1) .EQ. 0 ) THEN
                  NRVDEL = NRVDEL + 1
               END IF

            END DO

         END IF

C
C        Compute the base address of the current segment vector.
C
         SVBASE  =  JRSBAS + JSVBAS + (I-1)*SVSIZE

         IF (  ( RC .EQ. 0 ) .OR. ( NRVDEL .EQ. RC )  ) THEN
C
C           We're going to delete the current segment vector.  We'll
C           just skip over it without advancing our target pointers.
C
            NSVDEL   =  NSVDEL + 1

         ELSE IF ( NSVDEL .GT. 0 ) THEN
C
C           We need to shift the current segment vector to its
C           destination.
C
            CALL ZZEKSRD  ( SVBASE+1, SVBASE+SVSIZE, SEGVEC )
            CALL ZZEKSUPD ( VTARG+1,  VTARG+SVSIZE,  SEGVEC )

            VTARG = VTARG + SVSIZE

         ELSE
C
C           No segment vectors have been deleted yet.  We still must
C           update the target in case we shift vectors later on in this
C           loop.
C
            VTARG = VTARG + SVSIZE

         END IF

      END DO

C
C     At this point, we've compressed out the null segment vectors.
C     The next step is to compress out the row vector counts and row
C     vector pointers that corresponded to those segment vectors.  We
C     also want to remove the gap between the segment vectors and the
C     row vector pointer/count pairs.
C
C     We need to do this only if we deleted some segment vectors.
C
      IF ( NSVDEL .GT. 0 ) THEN

         NEWNSV =  NSV    - NSVDEL
         PTARG  =  JRSBAS + JSVBAS + (NEWNSV)*SVSIZE

         DO I = 1, NSV
C
C           The row count is RC.
C
            SVSIZE  =  NTAB
            CNTLOC  =  JRSBAS + JSVBAS + NSV*SVSIZE + 2*(I-1) + 2

            CALL ZZEKSRD ( CNTLOC, CNTLOC, RC )

            PTBASE  =  CNTLOC - 2

            IF ( RC .GT. 0 ) THEN
C
C              Shift the current row vector pointer and row vector
C              count.
C
               CALL ZZEKSRD  ( PTBASE+1, PTBASE+2, PCPAIR )
               CALL ZZEKSUPD ( PTARG+1,  PTARG+2,  PCPAIR )

               PTARG =  PTARG + 2

            END IF

         END DO

      ELSE

         NEWNSV = NSV

      END IF

C
C     Update the segment vector count.
C
      CALL ZZEKSUPD ( NSVLOC, NSVLOC, NEWNSV )

C
C     Remove any gaps that may exist between any of the row vectors,
C     or between the end of the segment vector's row vector counts
C     and base addresses and the first row vector.
C
C     The initial target location is the first element following the
C     last segment vector's row vector count.  RTARG is used as a base
C     address; it precedes this location by 1.
C
C     If we deleted any segment vectors, the segment vector pointers
C     embedded in the row vectors must change.  Make these updates
C     if necessary.
C
C
      RTARG   =  JRSBAS + JSVBAS + NEWNSV*(SVSIZE+2)

      DO I = 1, NEWNSV
C
C        Find the row count and row pointer for the current segment
C        vector.
C
         CNTLOC  =  JRSBAS + JSVBAS + NEWNSV*SVSIZE + 2*(I-1) + 2

         CALL ZZEKSRD ( CNTLOC, CNTLOC, RC )

         PTRLOC  =  CNTLOC - 1

C
C        Get the row vector set base pointer.  After capturing the
C        current value, we'll update this pointer to account for
C        the shifting of row vectors.
C
         CALL ZZEKSRD  ( PTRLOC, PTRLOC, SETBAS )

         RBASE   =  JRSBAS + SETBAS
         DELTA   =  RTARG  - RBASE

         CALL ZZEKSUPD ( PTRLOC, PTRLOC, SETBAS+DELTA )

C
C        Shift the row vectors for the current segment vector,
C        leaving behind the row vectors marked for deletion.
C
         NRVDEL = 0

         DO J = 1, RC

            CALL ZZEKSRD  ( RBASE+1, RBASE+RVSIZE, ROWVEC )

            IF ( ROWVEC(1) .EQ. 0 ) THEN
C
C              This row vector is to be deleted; don't copy it.
C
               RBASE  = RBASE  + RVSIZE
               NRVDEL = NRVDEL + 1

            ELSE
C
C              The segment vector pointer is base-relative.
C
               ROWVEC( RVSIZE )  =  JSVBAS + (I-1)*SVSIZE

               CALL ZZEKSUPD ( RTARG+1, RTARG+RVSIZE, ROWVEC )

               RBASE  =  RBASE + RVSIZE
               RTARG  =  RTARG + RVSIZE

            END IF

         END DO

C
C        Update the row count for the current segment vector, if
C        necessary.  Note that no segment vector will become empty
C        as a result of the row vector deletions we've done; we
C        already eliminated any segment vectors for which that
C        could happen, before we entered this loop.
C
         IF ( NRVDEL .GT. 0 ) THEN
            CALL ZZEKSUPD (  CNTLOC,  CNTLOC,  RC - NRVDEL  )
         END IF

      END DO

C
C     Update the total row count and size of the join row set.
C
      NR  =  0

      DO I = 1, NEWNSV

         CNTLOC  =  JRSBAS + JSVBAS + NEWNSV*SVSIZE + (I-1)*2 + 2

         CALL ZZEKSRD ( CNTLOC, CNTLOC, RC )

         NR      =  NR     + RC

      END DO

      NRLOC  =  JRSBAS   +   JRCIDX
      SIZE   =  JSVBAS   +   NEWNSV * ( SVSIZE + 2 )  +  NR * RVSIZE

      CALL ZZEKSUPD ( NRLOC,  NRLOC,  NR   )
      CALL ZZEKSUPD ( SIZLOC, SIZLOC, SIZE )

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!  Description: How to copy a namespace from a message to another.
!
!
program copy_namespace
   use eccodes
   implicit none
   integer  :: file1, file2, file3
   integer  :: igrib1, igrib2, igrib3

   call codes_open_file(file1, '../../data/reduced_latlon_surface.grib2', 'r')
   call codes_open_file(file2, '../../data/regular_latlon_surface.grib1', 'r')
   call codes_open_file(file3, 'out.grib_copy_namespace.grib', 'w')

   call codes_grib_new_from_file(file1, igrib1)
   call codes_grib_new_from_file(file2, igrib2)

   call codes_clone(igrib2, igrib3)

   call codes_copy_namespace(igrib1, 'geography', igrib3)

   call codes_write(igrib3, file3)

   call codes_close_file(file1)
   call codes_close_file(file2)
   call codes_close_file(file3)

end program copy_namespace

max_stars_repo_name: dvuckovic/eccodes

max_issues_repo_name: dvuckovic/eccodes

max_forks_repo_name: dvuckovic/eccodes


text:
      subroutine ksdfp (x,side,act,dfun,ndv,s,slope,y,p,h,hess,
     1                  isdflg)
      implicit double precision (a-h,o-z)
      dimension x(*),side(*),act(*),dfun(*),s(*),y(*),p(*),h(*)
      dimension hess(*)
c
c          routine to compute search direction using the
c          davidon-fletcher-powell method
c
c          author   - Gregory A. Wrenn
c          location - Lockheed Engineering and Sciences Co.
c                     144 Research Drive
c                     Hampton, Va. 23666
c
c          last modification - 19 July 1996
c
      if (isdflg .gt. 0) go to 40
c
c          reset approximate hessian matrix to initial condition
c
   10 continue
c
c          set h matrix to hess matrix
c
      isdflg = 0
      nh = ndv * (ndv + 1) / 2
      do 20 i = 1,nh
        h(i) = hess(i)
   20 continue
      do 30 i = 1,ndv
        act(i) = side(i)
   30 continue
      go to 80
c
c          dfp algorithm
c
   40 continue
c
c          check for validity of hessian matrix
c          in the case of side constraints
c
      do 50 i = 1,ndv
        if (side(i) .ne. act(i)) go to 10
   50 continue
c
c          compute vectors p and y
c
      do 60 i = 1,ndv
        y(i) = dfun(i) - y(i)
        p(i) = x(i) - p(i)
   60 continue
c
c          compute update to hessian matrix
c
      call kshmul (h,y,s,ndv)
      call ksvprd (p,y,sigma,ndv)
      call ksvprd (y,s,tau,ndv)
      if (abs(sigma) .lt. 1.0e-08) sigma = 1.0e-08
      if (abs(tau) .lt. 1.0e-08) tau = 1.0e-08
      sigma = 1.0 / sigma
      tau = 1.0 / tau
c
      k = 0
      do 70 i = 1,ndv
      do 70 j = 1,i
        k = k + 1
        h(k) = h(k) + sigma * p(i) * p(j) - tau * s(i) * s(j)
   70 continue
c
   80 continue
      call kshmul (h,dfun,s,ndv)
c
c          negate s-vector
c
      smax = 1.0e-08
      do 90 i = 1,ndv
        s(i) = -s(i)
        si = abs(s(i))
        if (si .gt. smax) smax = si
   90 continue
c
c          normalize search direction to unit maximum
c          and save previous gradients in y-vector
c          and previous design variables in p-vector
c
      do 100 i = 1,ndv
        s(i) = s(i) / smax
        y(i) = dfun(i)
        p(i) = x(i)
  100 continue
c
c          compute slope
c
      call ksvprd (s,dfun,slope,ndv)
      if (slope .gt. 0.0 .and. isdflg .gt. 0) go to 10
      if (slope .gt. 0.0) slope = 0.0
c
      return
      end

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
C**********************************************************************
C
C     Copyright (C) 1992  Roland W. Freund and Noel M. Nachtigal
C     All rights reserved.
C
C     This code is part of a copyrighted package.  For details, see the
C     file `cpyrit.doc' in the top-level directory.
C
C     *****************************************************************
C     ANY USE OF  THIS CODE CONSTITUTES ACCEPTANCE OF  THE TERMS OF THE
C                             COPYRIGHT NOTICE
C     *****************************************************************
C
C**********************************************************************
C
C     This  file  contains  the  routine  for  the  QMR  algorithm  for
C     unsymmetric matrices, using the  three-term recurrence variant of
C     the Lanczos algorithm without look-ahead.
C
C**********************************************************************
C
      SUBROUTINE ZUQMX (NDIM,NLEN,NLIM,VECS,TOL,INFO)
C
C     Purpose:
C     This subroutine uses  the QMR algorithm to solve  linear systems.
C     It runs  the algorithm  to convergence or until  a user-specified
C     limit on the number of  iterations is reached.
C
C     The  code is  set up  to solve  the system  A x = b  with initial
C     guess x_0 = 0.  Here  A x = b  denotes the preconditioned system,
C     and  it  is  connected with the  original system as follows.  Let
C     B y = c be the original unpreconditioned system to be solved, and
C     let y_0 be an arbitrary initial guess for its solution.  Then:
C          A x = b, where  A = M_1^{-1} B M_2^{-1},
C                          x = M_2 (y - y_0), b = M_1^{-1} (c - B y_0).
C     Here M = M_1 M_2 is the preconditioner.
C
C     To recover the final iterate y_n for  the original system B y = c
C     from the final iterate x_n for the preconditioned system A x = b,
C     set
C               y_n = y_0 + M_2^{-1} x_n.
C
C     The implementation does not have look-ahead, so it is less robust
C     than the full version.
C
C     Parameters:
C     For a description of  the parameters, see the file `zuqmx.doc' in
C     the current directory.
C
C     External routines used:
C     double precision dlamch(ch)
C        LAPACK routine, computes machine-related constants.
C     double precision dznrm2(n,x,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine zaxpy(n,a,x,incx,y,incy)
C        BLAS-1 routine, computes y = a * x + y.
C     subroutine zaxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     double precision zdotu(n,x,incx,y,incy)
C        BLAS-1 routine, computes y' * x.
C     double precision zqmxom(n)
C        User-supplied routine, specifies the QMR scaling factors.
C     subroutine zrandn(n,x,seed)
C        Library routine, fills x with random numbers.
C     subroutine zrotg(a,b,cos,sin)
C        BLAS-1 routine, computes the Givens rotation which rotates the
C        vector [a; b] into [ sqrt(a**2 + b**2); 0 ].
C
C     Noel M. Nachtigal
C     May 25, 1993
C
C**********************************************************************
C
      INTRINSIC CDABS, DABS, DBLE, DCMPLX, DCONJG, DMAX1, DSQRT, MAX0
      INTRINSIC MOD
c
      EXTERNAL DLAMCH, DZNRM2, ZAXPBY, ZDOTU, ZRANDN, ZROTG, ZUQMXO
C
      integer    INFO(4)
      integer    NDIM
      integer    NLEN
      integer    NLIM
c
      real*8     DLAMCH
      real*8     DZNRM2
      real*8     ZUQMXO
      real*8     TOL
c
      complex*16 ZDOTU
      complex*16 VECS(NDIM,9)
C
C     Miscellaneous parameters.
C
      complex*16  ZONE
      complex*16  ZZERO
      PARAMETER (ZONE = (1.0D0,0.0D0),ZZERO = (0.0D0,0.0D0))
c
      real*8      DHUN
      real*8      DONE
      real*8      DTEN
      real*8      DZERO
      PARAMETER (DHUN = 1.0D2,DONE = 1.0D0,DTEN = 1.0D1,DZERO = 0.0D0)
C
C     Local variables, permanent.
C
      INTEGER IERR, ISN, ISNM1, ISNM2, IVN, IVNM1, IVNP1, IWN, IWNM1
      SAVE    IERR, ISN, ISNM1, ISNM2, IVN, IVNM1, IVNP1, IWN, IWNM1
c
      INTEGER IWNP1, N, RETLBL, TF, TRES, VF
      SAVE    IWNP1, N, RETLBL, TF, TRES, VF

      complex*16     DN, DNM1, DNP1, RHSN, SCSN, SCSNM1, SINN, SINNM1
      SAVE           DN, DNM1, DNP1, RHSN, SCSN, SCSNM1, SINN, SINNM1

      real*8           COSN, COSNM1, CSIN, CSINP1, MAXOMG, OMGN
      SAVE             COSN, COSNM1, CSIN, CSINP1, MAXOMG, OMGN

      real*8           OMGNP1, R0, RESN, RHON, RHONP1, SCVN, SCVNP1
      SAVE             OMGNP1, R0, RESN, RHON, RHONP1, SCVN, SCVNP1

      real*8           SCWN, SCWNP1, TMAX, TMIN, TNRM, UCHK, UNRM
      SAVE             SCWN, SCWNP1, TMAX, TMIN, TNRM, UCHK, UNRM
C
C     Local variables, transient.
C
      integer       INIT
      integer       REVCOM
c
      complex*16    RHN
      complex*16    RHNM1
      complex*16    RHNM2
      complex*16    RHNP1
      complex*16    RHSNP1
      complex*16    ZTMP
c
      real*8        DTMP
      real*8        SCVNM1
      real*8        SCWNM1
C
C     Initialize some of the permanent variables.
C
      DATA RETLBL /0/
C
C     Check the reverse communication flag to see where to branch.
C        REVCOM   RETLBL      Comment
C           0        0    first call, go to label 10
C           1       30    returning from AXB, go to label 30
C           1       60    returning from AXB, go to label 60
C           2       40    returning from ATXB, go to label 40
C
      REVCOM  = INFO(2)
      INFO(2) = 0
      IF (REVCOM.EQ.0) THEN
         N = 0
         IF (RETLBL.EQ.0) GO TO 10
      ELSE IF (REVCOM.EQ.1) THEN
         IF (RETLBL.EQ.30) THEN
            GO TO 30
         ELSE IF (RETLBL.EQ.60) THEN
            GO TO 60
         END IF
      ELSE IF (REVCOM.EQ.2) THEN
         IF (RETLBL.EQ.40) GO TO 40
      END IF
      IERR = 1
      GO TO 80
C
C     Check whether the inputs are valid.
C
 10   IERR = 0
      IF (NDIM.LT.1)        IERR = 2
      IF (NLEN.LT.1)        IERR = 2
      IF (NLIM.LT.1)        IERR = 2
      IF (NLEN.GT.NDIM)     IERR = 2
      IF (IERR.NE.0) GO TO 80
C
C     Extract from INFO the output units TF and VF, the true residual
C     flag TRES, and the left starting vector flag INIT.
C
      VF   = MAX0(INFO(1),0)
      INIT = VF / 100000
      VF   = VF - INIT * 100000
      TRES = VF / 10000
      VF   = VF - TRES * 10000
      TF   = VF / 100
      VF   = VF - TF * 100
C
C     Extract and check the various tolerances.
C
      TNRM = DLAMCH('E') * DTEN
      TMIN = DSQRT(DSQRT(DLAMCH('S')))
      TMAX = DONE / TMIN
      IF (TOL.LE.DZERO) TOL = DSQRT(DLAMCH('E'))
C
C     Start the trace messages and convergence history.
C
      IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') 0, DONE, DONE
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') 0, DONE, DONE
C
C     Initialize the wrapped indices.
C
      ISNM1 = 5
      ISN   = ISNM1
      IVN   = 3
      IVNP1 = IVN
      IWN   = 8
      IWNP1 = IWN
C
C     Set x_0 = 0 and compute the norm of the initial residual.
C
c     CALL ZAXPBY (NLEN,VECS(1,IVN),ZONE,VECS(1,2),ZZERO,VECS(1,IVN))
c     CALL ZAXPBY (NLEN,VECS(1,1),ZZERO,VECS(1,1),ZZERO,VECS(1,1))
      CALL Zcopy (NLEN,VECS(1,2),1,VECS(1,IVN),1)
      CALL Zscal (NLEN,zzero,VECS(1,1),1)
      R0 = DZNRM2(NLEN,VECS(1,IVN),1)
      IF ((TOL.GE.DONE).OR.(R0.EQ.DZERO)) GO TO 80
C
C     Check whether the auxiliary vector must be supplied.
C
      IF (INIT.EQ.0) CALL ZRANDN (NLEN,VECS(1,IWN),1)
C
C     Compute scale factors and check for invariant subspaces.
C
      SCVNP1 = R0
      SCWNP1 = DZNRM2(NLEN,VECS(1,IWN),1)
      IF (SCVNP1.LT.TNRM) IERR = IERR + 16
      IF (SCWNP1.LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GO TO 80
      DNP1 = ZDOTU(NLEN,VECS(1,IWN),1,VECS(1,IVN),1) / ( SCVNP1 *
     $SCWNP1 )
      IF ((SCVNP1.GE.TMAX).OR.(SCVNP1.LE.TMIN)) THEN
         ZTMP = DCMPLX(DONE / SCVNP1,DZERO)
c        CALL ZAXPBY (NLEN,VECS(1,IVN),ZTMP,VECS(1,IVN),ZZERO,VECS(1,
c    $IVN))
         CALL Zscal (NLEN,ZTMP,VECS(1,IVN),1)
         SCVNP1 = DONE
      END IF
      IF ((SCWNP1.GE.TMAX).OR.(SCWNP1.LE.TMIN)) THEN
         ZTMP = DCMPLX(DONE / SCWNP1,DZERO)
c        CALL ZAXPBY (NLEN,VECS(1,IWN),ZTMP,VECS(1,IWN),ZZERO,VECS(1,
c    $IWN))
         CALL Zscal (NLEN,ZTMP,VECS(1,IWN),1)
         SCWNP1 = DONE
      END IF
      RHONP1 = SCVNP1
      CSINP1 = SCWNP1
      SCVNP1 = DONE / SCVNP1
      SCWNP1 = DONE / SCWNP1
C
C     Initialize the variables.
C
      N      = 1
      DN     = ZONE
      COSN   = DONE
      RESN   = DONE
      COSNM1 = DONE
      OMGN   = DZERO
      SCSN   = DZERO
      SCVN   = DZERO
      SCWN   = DZERO
      SCSNM1 = DZERO
      SINN   = ZZERO
      SINNM1 = ZZERO
      OMGNP1 = ZUQMXO(N)
      RHSN   = OMGNP1 * R0
      MAXOMG = DONE / OMGNP1
C
C     This is one step of the classical Lanczos algorithm.
C
 20   IVNM1 = IVN
      IVN   = IVNP1
      IVNP1 = MOD(N,2) + 3
      IWNM1 = IWN
      IWN   = IWNP1
      IWNP1 = MOD(N,2) + 8
C
C     Check whether D_n is nonsingular.
C
      DNM1 = DN
      DN   = DNP1
      IF (CDABS(DN).EQ.DZERO) THEN
         IERR = 8
         GO TO 80
      END IF
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,IVN),VECS(1,7))
C
      INFO(2) = 1
      INFO(3) = IVN
      INFO(4) = 7
      RETLBL  = 30
      RETURN
 30   RETLBL = 0
C
C     Compute H_{n-1,n} and build part of the vector v_{n+1}.
C
      SCVNM1 = SCVN
      CSIN   = CSINP1
      SCVN   = SCVNP1
      ZTMP   = CSIN * DN / DNM1 * SCVNM1 / SCVN
      CALL ZAXPBY (NLEN,VECS(1,IVNP1),ZONE,VECS(1,7),-ZTMP,VECS(1,IVNM1)
     $)
C
C     Have the caller carry out ATXB, then return here.
C        CALL ATXB (VECS(1,IWN),VECS(1,7))
C
      INFO(2) = 2
      INFO(3) = IWN
      INFO(4) = 7
      RETLBL  = 40
      RETURN
 40   RETLBL = 0
C
C     Build part of the vector w_{n+1}.
C
      SCWNM1 = SCWN
      RHON   = RHONP1
      SCWN   = SCWNP1
      ZTMP   = RHON * DN / DNM1 * SCWNM1 / SCWN
      CALL ZAXPBY (NLEN,VECS(1,IWNP1),ZONE,VECS(1,7),-ZTMP,VECS(1,IWNM1)
     $)
C
C     Compute H_{nn} and finish the new vectors.
C
      RHN = SCVN * SCWN * ZDOTU(NLEN,VECS(1,IWN),1,VECS(1,IVNP1),1) / DN
c     CALL ZAXPBY (NLEN,VECS(1,IVNP1),ZONE,VECS(1,IVNP1),-RHN,VECS(1,
c    $IVN))
c     CALL ZAXPBY (NLEN,VECS(1,IWNP1),ZONE,VECS(1,IWNP1),-RHN,VECS(1,
c    $IWN))
      CALL ZAXPY (NLEN,-RHN,VECS(1,IVN),1,VECS(1,IVNP1),1)
      CALL ZAXPY (NLEN,-RHN,VECS(1,IWN),1,VECS(1,IWNP1),1)
C
C     Compute scale factors and check for invariant subspaces.
C
      IERR   = 0
      SCVNP1 = DZNRM2(NLEN,VECS(1,IVNP1),1)
      SCWNP1 = DZNRM2(NLEN,VECS(1,IWNP1),1)
      RHONP1 = SCVN * SCVNP1
      CSINP1 = SCWN * SCWNP1
      RHNP1  = RHONP1
      IF (SCVNP1.LT.TNRM) IERR = IERR + 16
      IF (SCWNP1.LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GO TO 50
      DNP1   = ZDOTU(NLEN,VECS(1,IWNP1),1,VECS(1,IVNP1),1) / ( SCVNP1 *
     $SCWNP1 )
      IF ((SCVNP1.GE.TMAX).OR.(SCVNP1.LE.TMIN)) THEN
         ZTMP = DCMPLX(DONE / SCVNP1,DZERO)
c        CALL ZAXPBY (NLEN,VECS(1,IVNP1),ZTMP,VECS(1,IVNP1),ZZERO,
c    $VECS(1,IVNP1))
         CALL Zscal (NLEN,ZTMP,VECS(1,IVNP1),1)
         SCVNP1 = DONE
      END IF
      IF ((SCWNP1.GE.TMAX).OR.(SCWNP1.LE.TMIN)) THEN
         ZTMP = DCMPLX(DONE / SCWNP1,DZERO)
c        CALL ZAXPBY (NLEN,VECS(1,IWNP1),ZTMP,VECS(1,IWNP1),ZZERO,
c    $VECS(1,IWNP1))
         CALL Zscal (NLEN,ZTMP,VECS(1,IWNP1),1)
         SCWNP1 = DONE
      END IF
      SCVNP1 = DONE / SCVNP1
      SCWNP1 = DONE / SCWNP1
C
C     The QMR code starts here.
C     Multiply the new column by the previous omega's.
C     Get the next scaling factor omega(i) and update MAXOMG.
C
 50   RHNM1  = CSIN * DN * OMGN / DNM1
      OMGN   = OMGNP1
      RHN    = OMGN * RHN
      OMGNP1 = ZUQMXO(N+1)
      RHNP1  = OMGNP1 * RHNP1
      MAXOMG = DMAX1(MAXOMG,DONE/OMGN)
C
C     Apply the previous rotations.
C
      RHNM2  = SINNM1 * RHNM1
      RHNM1  = COSNM1 * RHNM1
      COSNM1 = COSN
      SINNM1 = SINN
      ZTMP   = RHNM1
      RHNM1  =  COSNM1 * ZTMP + SINNM1 * RHN
      RHN    = -DCONJG(SINNM1) * ZTMP + COSNM1 * RHN
C
C     Compute the rotation for the last element (this also applies it).
C
      CALL ZROTG (RHN,RHNP1,COSN,SINN)
C
C     Apply the new rotation to the right-hand side vector.
C
      RHSNP1 = -DCONJG(SINN) * RHSN
      RHSN   =  COSN * RHSN
C
C     Compute the next search direction s_n.
C
      ISNM2  = ISNM1
      ISNM1  = ISN
      ISN    = MOD(N-1,2) + 5
      ZTMP   = SCSNM1 * RHNM2 / SCVN
      CALL ZAXPBY (NLEN,VECS(1,ISN),ZONE,VECS(1,IVN),-ZTMP,VECS(1,ISNM2)
     $)
      SCSNM1 = SCSN
      ZTMP   = SCSNM1 * RHNM1 / SCVN
c     CALL ZAXPBY (NLEN,VECS(1,ISN),ZONE,VECS(1,ISN),-ZTMP,VECS(1,ISNM1)
c    $)
      CALL ZAXPY (NLEN,-ZTMP,VECS(1,ISNM1),1,VECS(1,ISN),1)
      SCSN   = SCVN / RHN
C
C     Compute the new QMR iterate, then scale the search direction.
C
      ZTMP = SCSN * RHSN
c     CALL ZAXPBY (NLEN,VECS(1,1),ZONE,VECS(1,1),ZTMP,VECS(1,ISN))
      CALL ZAXPY (NLEN,ZTMP,VECS(1,ISN),1,VECS(1,1),1)
      DTMP = CDABS(SCSN)
      IF ((DTMP.GE.TMAX).OR.(DTMP.LE.TMIN)) THEN
c        CALL ZAXPBY (NLEN,VECS(1,ISN),SCSN,VECS(1,ISN),ZZERO,VECS(1,
c    $ISN))
         CALL Zscal (NLEN,SCSN,VECS(1,ISN),1)
         SCSN = ZONE
      END IF
C
C     Compute the residual norm upper bound.
C     If the scaled upper bound is within one order of magnitude of the
C     target convergence norm, compute the true residual norm.
C
      RHSN = RHSNP1
      UNRM = DSQRT(DBLE(N+1)) * MAXOMG * CDABS(RHSNP1) / R0
      UCHK = UNRM
      IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.DTEN).AND.(N.LT.NLIM)) GO TO 70
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,1),VECS(1,7))
C
      INFO(2) = 1
      INFO(3) = 1
      INFO(4) = 7
      RETLBL  = 60
      RETURN
 60   RETLBL = 0
      CALL ZAXPBY (NLEN,VECS(1,7),ZONE,VECS(1,2),-ZONE,VECS(1,7))
      RESN = DZNRM2(NLEN,VECS(1,7),1) / R0
      UCHK = RESN
C
C     Output the trace messages and convergence history.
C
 70   IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') N, UNRM, RESN
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') N, UNRM, RESN
C
C     Check for convergence or termination.  Stop if:
C         1. algorithm converged;
C         2. there is an error condition;
C         3. the residual norm upper bound is smaller than the computed
C     residual norm by a factor of at least 100;
C         4. algorithm exceeded the iterations limit.
C
      IF (RESN.LE.TOL) THEN
         IERR = 0
         GO TO 80
      ELSE IF (IERR.NE.0) THEN
         GO TO 80
      ELSE IF (UNRM.LT.UCHK/DHUN) THEN
         IERR = 4
         GO TO 80
      ELSE IF (N.GE.NLIM) THEN
         IERR = 4
         GO TO 80
      END IF
C
C     Update the running counter.
C
      N = N + 1
      GO TO 20
C
C     That's all.
C
 80   ZLIM = N
      RETLBL  = 0
      INFO(1) = IERR
C
      RETURN
      END
C
C**********************************************************************
C
      FUNCTION ZUQMXO (I)
C
C     Purpose:
C     Returns the scaling parameter OMEGA(I).
C
C     Parameters:
C     I = the index of the parameter OMEGA (input).
C
C     Noel M. Nachtigal
C     March 30, 1993
C
C**********************************************************************
C
      integer  I
c
      real*8   ZUQMXO
C
      ZUQMXO = 1.0D0
C
      RETURN
      END
C
C**********************************************************************

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: hkershaw-brown/MuST


text:
      SUBROUTINE JULCAL(MONTH,IDAY,IYEAR,JDAY)
!     JULCAL CALCULATES AND RETURNS THE CALENDAR DATE GIVEN THE JULIAN DAY AND
!  YEAR.  LEAP YEARS ARE ACCOUNTED FOR.
!
!  ARGUMENTS:
!    MONTH  - THE MONTH OF THE YEAR CALCULATED (AND RETURNED) BY JULCAL.
!    IDAY   - THE DAY OF THE MONTH CALCULATED (AND RETURNED) BY JULCAL.
!    IYEAR  - THE YEAR OF THE JULIAN DAY TO BE CONVERTED.
!    JDAY   - THE JULIAN DAY TO BE CONVERTED.
!
!  CONVERTED FROM THE IBM1800 SYSTEM BY PAUL HENKART, SCRIPPS INSTITUTION OF
!  OCEANOGRAPHY, 28 OCTOBER 1980
!
      DIMENSION MTHDY(12)
      DATA MTHDY/31,28,31,30,31,30,31,31,30,31,30,31/
!
      ISUM=0
!
!     CHECK FOR LEAP YEAR
!
!      IF(IYEAR-IYEAR/4*4)2,1,2
      IF(IYEAR-IYEAR/4*4 .EQ. 0 ) MTHDY(2)=29
!
!     SUM DAYS/ MONTH
!
2     DO 3 I=1,12
      ISUM=ISUM+MTHDY(I)
      ii = i
!
!     COMPARE AGAINST JULIAN DAY
!
!      IF(ISUM-JDAY) 3,4,4
      IF(ISUM-JDAY .GE. 0 ) GOTO 4
3     CONTINUE
!
!     JULIAN DAY FALLS WITHIN MONTH I, SO SUBTRACT OUT MONTH I FOR DATE
!
4     MONTH=II
      IDAY=JDAY+MTHDY(II)-ISUM
!
!     RESET LEAP YEAR COUNT
!
      MTHDY(2)=28
      RETURN
      END


! ============================================================================
! Name        : mo_eq_lv.f90
! Author      : Giovanni Dalmasso
! Version     : 1.0
! Copyright   : Department of Computational Hydrosystems CHS - UZF Leipzig
! Description : Lotka Volterra equations
! ============================================================================

module mo_eq_lv

    ! ------------------------------------------------------------------
    ! The Lotka Volterra equations, also known as the predator-prey equations, are a pair of first-order, non-linear,
    ! differential equations frequently used to describe the dynamics of biological systems in which two species interact,
    ! one a predator and one its prey. They evolve in time according to the following of equations
    ! ------------------------------------------------------------------

    implicit none

    private

    public :: LV_eqn, LV_eqn_sp, LV_eqn_dp, &
         LV_eqn_para, LV_eqn_para_sp, LV_eqn_para_dp, &
         deriv_testRB_sp, jacobn_testRB_sp, &
         deriv_testRB_dp, jacobn_testRB_dp, &
         deriv2_testRB_sp, jacobn2_testRB_sp, &
         deriv2_testRB_dp, jacobn2_testRB_dp, &
         deriv2_testRB_para_sp, jacobn2_testRB_para_sp, &
         deriv2_testRB_para_dp, jacobn2_testRB_para_dp

    ! Interfaces for single and double precision routines
    interface LV_eqn
        module procedure LV_eqn_sp, LV_eqn_dp
    end interface

    interface LV_eqn_para
        module procedure LV_eqn_para_sp, LV_eqn_para_dp
    end interface

contains

    ! SINGLE PRECISION
    subroutine LV_eqn_sp( x, y, dydx )    ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                       ! time
        real(sp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(sp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(sp), parameter :: epsilonn = 0.5_sp        ! regulates the relationship between the two species
                                                        ! keeping the optimum balance

        real(sp), parameter :: mu = 2._sp               ! represents the ability to reproduce of the first species
                                                        ! with respect to the destructive capacity of the second one

        real(sp), parameter :: eta = -1._sp             ! regulates the relationship between the influences of nutrition
                                                        ! related reproductive capacity

        ! ============================================================================
        ! System of ODE
        !
        ! y(1) --> predators
        ! y(2) --> preys
        ! does not depend on time  --> 0.0 * x (to avoid warnings of unused variables)
        ! ============================================================================

        dydx(1) = y(1) - epsilonn*mu*y(2)        + 0.0_sp*x
        dydx(2) = (mu/epsilonn)*y(1) + eta*y(2)  + 0.0_sp*x

    end subroutine LV_eqn_sp

    ! DOUBLE PRECISION
    subroutine LV_eqn_dp( x, y, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                       ! time
        real(dp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(dp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(dp), parameter :: epsilonn = 0.5_dp        ! regulates the relationship between the two species
                                                        ! keeping the optimum balance

        real(dp), parameter :: mu = 2._dp               ! represents the ability to reproduce of the first species
                                                        ! with respect to the destructive capacity of the second one

        real(dp), parameter :: eta = -1._dp             ! regulates the relationship between the influences of nutrition
                                                        ! related reproductive capacity

        ! ============================================================================
        ! System of ODE
        !
        ! y(1) --> predators
        ! y(2) --> preys
        ! does not depend on time  --> 0.0 * x (to avoid warnings of unused variables)
        ! ============================================================================

        dydx(1) = y(1) - epsilonn*mu*y(2)        + 0.0_dp*x
        dydx(2) = (mu/epsilonn)*y(1) + eta*y(2)  + 0.0_dp*x

    end subroutine LV_eqn_dp

!-----------
! with parameter input
!-----------

    ! SINGLE PRECISION
    subroutine LV_eqn_para_sp( x, y, para, dydx )    ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                       ! time
        real(sp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(sp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(sp), dimension(:), intent(in) :: para      ! parameter

!        real(sp), parameter :: epsilonn = 0.5_sp        ! regulates the relationship between the two species
!                                                        ! keeping the optimum balance

!        real(sp), parameter :: mu = 2._sp               ! represents the ability to reproduce of the first species
!                                                        ! with respect to the destructive capacity of the second one

!        real(sp), parameter :: eta = -1._sp             ! regulates the relationship between the influences of nutrition
!                                                        ! related reproductive capacity

        ! ============================================================================
        ! System of ODE
        !
        ! y(1) --> predators
        ! y(2) --> preys
        ! does not depend on time  --> 0.0 * x (to avoid warnings of unused variables)
        ! ============================================================================

        dydx(1) = y(1) - para(1)*para(2)*y(2)           + 0.0_sp*x
        dydx(2) = (para(2)/para(1))*y(1) + para(3)*y(2) + 0.0_sp*x

    end subroutine LV_eqn_para_sp

    ! DOUBLE PRECISION
    subroutine LV_eqn_para_dp( x, y, para, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                       ! time
        real(dp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(dp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(dp), dimension(:), intent(in) :: para      ! parameter

!        real(dp), parameter :: epsilonn = 0.5_dp        ! regulates the relationship between the two species
!                                                        ! keeping the optimum balance

!        real(dp), parameter :: mu = 2._dp               ! represents the ability to reproduce of the first species
!                                                        ! with respect to the destructive capacity of the second one

!        real(dp), parameter :: eta = -1._dp             ! regulates the relationship between the influences of nutrition
!                                                        ! related reproductive capacity

        ! ============================================================================
        ! System of ODE
        !
        ! y(1) --> predators
        ! y(2) --> preys
        ! does not depend on time  --> 0.0 * x (to avoid warnings of unused variables)
        ! ============================================================================

        dydx(1) = y(1) - para(1)*para(2)*y(2)           + 0.0_dp*x
        dydx(2) = (para(2)/para(1))*y(1) + para(3)*y(2) + 0.0_dp*x

    end subroutine LV_eqn_para_dp

!-----------
! stiff sets of ODEs
!-----------

!SINGLE PRECISION
    subroutine deriv_testRB_sp( x, y, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                       ! time
        real(sp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(sp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(sp) :: xtmp

        xtmp = x

        dydx(1) = 998.0_sp*y(1) + 1998.0_sp*y(2)
        dydx(2) = -999.0_sp*y(1)- 1999.0_sp*y(2)

    end subroutine deriv_testRB_sp

    subroutine jacobn_testRB_sp( x, y, dfdx, dfdy )   ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                           ! time
        real(sp), dimension(:),     intent(in)  :: y        ! unknowns of the equations
        real(sp), dimension(:),     intent(out) :: dfdx     ! derivatives of f for x
        real(sp), dimension(:,:),   intent(out) :: dfdy     ! derivatives of f for y

        real(sp) :: xtmp
        real(sp), dimension(size(y,1)) :: ytmp

        xtmp = x
        ytmp = y

        dfdx(1)     = 0.0_sp
        dfdx(2)     = 0.0_sp

        dfdy(1,1)   = 998.0_sp
        dfdy(1,2)   = 1998.0_sp
        dfdy(2,1)   = -999.0_sp
        dfdy(2,2)   = 1999.0_sp

    end subroutine jacobn_testRB_sp

    subroutine deriv2_testRB_sp( x, y, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                       ! time
        real(sp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(sp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(sp) :: xtmp

        xtmp = x

        dydx(1) = (-0.013_sp-1000.0_sp*y(3))*y(1)
        dydx(2) = -2500.0_sp*y(3)*y(2)
        dydx(3) = -0.013_sp*y(1) + (-1000.0_sp*y(1)-2500.0_sp*y(2))*y(3)

    end subroutine deriv2_testRB_sp

    subroutine jacobn2_testRB_sp( x, y, dfdx, dfdy )   ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                           ! time
        real(sp), dimension(:),     intent(in)  :: y        ! unknowns of the equations
        real(sp), dimension(:),     intent(out) :: dfdx     ! derivatives of f for x
        real(sp), dimension(:,:),   intent(out) :: dfdy     ! derivatives of f for y

        real(sp) :: xtmp
        real(sp), dimension(size(y,1)) :: ytmp

        xtmp = x
        ytmp = y

        dfdx(1)     = 0.0_sp
        dfdx(2)     = 0.0_sp
        dfdx(3)     = 0.0_sp

        dfdy(1,1)   = -0.013_sp-1000.0_sp*y(3)
        dfdy(1,2)   = 0.0_sp
        dfdy(1,3)   = -1000.0_sp*y(1)
        dfdy(2,1)   = 0.0_sp
        dfdy(2,2)   = -2500.0_sp*y(3)
        dfdy(2,3)   = -2500.0_sp*y(2)
        dfdy(3,1)   = -0.013_sp-1000.0_sp*y(3)
        dfdy(3,2)   = -2500.0_sp*y(3)
        dfdy(3,3)   = -1000.0_sp*y(1)-2500.0_sp*y(2)

    end subroutine jacobn2_testRB_sp

    subroutine deriv2_testRB_para_sp( x, y, para, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                       ! time
        real(sp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(sp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(sp), dimension(:), intent(in) :: para      ! parameter

        real(sp) :: xtmp

        xtmp = x

        dydx(1) = (-0.013_sp-para(1)*y(3))*y(1)
        dydx(2) = -2500.0_sp*y(3)*y(2)
        dydx(3) = -0.013_sp*y(1) + (-1000.0_sp*y(1)-2500.0_sp*y(2))*y(3)

    end subroutine deriv2_testRB_para_sp

    subroutine jacobn2_testRB_para_sp( x, y, para, dfdx, dfdy )   ! returns derivatives dydx at x

        use mo_kind, only : sp

        implicit none

        real(sp), intent(in) :: x                           ! time
        real(sp), dimension(:),     intent(in)  :: y        ! unknowns of the equations
        real(sp), dimension(:),     intent(out) :: dfdx     ! derivatives of f for x
        real(sp), dimension(:,:),   intent(out) :: dfdy     ! derivatives of f for y

        real(sp), dimension(:),     intent(in)  :: para      ! parameter

        real(sp) :: xtmp
        real(sp), dimension(size(y,1)) :: ytmp

        xtmp = x
        ytmp = y

        dfdx(1)     = 0.0_sp
        dfdx(2)     = 0.0_sp
        dfdx(3)     = 0.0_sp

        dfdy(1,1)   = -0.013_sp-para(1)*y(3)
        dfdy(1,2)   = 0.0_sp
        dfdy(1,3)   = -1000.0_sp*y(1)
        dfdy(2,1)   = 0.0_sp
        dfdy(2,2)   = -2500.0_sp*y(3)
        dfdy(2,3)   = -2500.0_sp*y(2)
        dfdy(3,1)   = -0.013_sp-1000.0_sp*y(3)
        dfdy(3,2)   = -2500.0_sp*y(3)
        dfdy(3,3)   = -1000.0_sp*y(1)-2500.0_sp*y(2)

    end subroutine jacobn2_testRB_para_sp

!DOUBLE PRECISION
    subroutine deriv_testRB_dp( x, y, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                       ! time
        real(dp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(dp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(dp) :: xtmp

        xtmp = x

        dydx(1) = 998.0_dp*y(1) + 1998.0_dp*y(2)
        dydx(2) = -999.0_dp*y(1)- 1999.0_dp*y(2)

    end subroutine deriv_testRB_dp

    subroutine jacobn_testRB_dp( x, y, dfdx, dfdy )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                           ! time
        real(dp), dimension(:),     intent(in)  :: y        ! unknowns of the equations
        real(dp), dimension(:),     intent(out) :: dfdx     ! derivatives of f for x
        real(dp), dimension(:,:),   intent(out) :: dfdy     ! derivatives of f for y

        real(dp) :: xtmp
        real(dp), dimension(size(y,1)) :: ytmp

        xtmp = x
        ytmp = y

        dfdx(1)     = 0.0_dp
        dfdx(2)     = 0.0_dp

        dfdy(1,1)   = 998.0_dp
        dfdy(1,2)   = 1998.0_dp
        dfdy(2,1)   = -999.0_dp
        dfdy(2,2)   = 1999.0_dp

    end subroutine jacobn_testRB_dp

    subroutine deriv2_testRB_dp( x, y, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                       ! time
        real(dp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(dp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(dp) :: xtmp

        xtmp = x

        dydx(1) = (-0.013_dp-1000.0_dp*y(3))*y(1)
        dydx(2) = -2500.0_dp*y(3)*y(2)
        dydx(3) = -0.013_dp*y(1) + (-1000.0_dp*y(1)-2500.0_dp*y(2))*y(3)

    end subroutine deriv2_testRB_dp

    subroutine jacobn2_testRB_dp( x, y, dfdx, dfdy )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                           ! time
        real(dp), dimension(:),     intent(in)  :: y        ! unknowns of the equations
        real(dp), dimension(:),     intent(out) :: dfdx     ! derivatives of f for x
        real(dp), dimension(:,:),   intent(out) :: dfdy     ! derivatives of f for y

        real(dp) :: xtmp
        real(dp), dimension(size(y,1)) :: ytmp

        xtmp = x
        ytmp = y

        dfdx(1)     = 0.0_dp
        dfdx(2)     = 0.0_dp
        dfdx(3)     = 0.0_dp

        dfdy(1,1)   = -0.013_dp-1000.0_dp*y(3)
        dfdy(1,2)   = 0.0_dp
        dfdy(1,3)   = -1000.0_dp*y(1)
        dfdy(2,1)   = 0.0_dp
        dfdy(2,2)   = -2500.0_dp*y(3)
        dfdy(2,3)   = -2500.0_dp*y(2)
        dfdy(3,1)   = -0.013_dp-1000.0_dp*y(3)
        dfdy(3,2)   = -2500.0_dp*y(3)
        dfdy(3,3)   = -1000.0_dp*y(1)-2500.0_dp*y(2)

    end subroutine jacobn2_testRB_dp

    subroutine deriv2_testRB_para_dp( x, y, para, dydx )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                       ! time
        real(dp), dimension(:), intent(in) :: y         ! unknowns of the equations
        real(dp), dimension(:), intent(out) :: dydx     ! derivatives of y

        real(dp), dimension(:), intent(in) :: para      ! parameter

        real(dp) :: xtmp

        xtmp = x

        dydx(1) = (-0.013_dp-para(1)*y(3))*y(1)
        dydx(2) = -2500.0_dp*y(3)*y(2)
        dydx(3) = -0.013_dp*y(1) + (-1000.0_dp*y(1)-2500.0_dp*y(2))*y(3)

    end subroutine deriv2_testRB_para_dp

    subroutine jacobn2_testRB_para_dp( x, y, para, dfdx, dfdy )   ! returns derivatives dydx at x

        use mo_kind, only : dp

        implicit none

        real(dp), intent(in) :: x                           ! time
        real(dp), dimension(:),     intent(in)  :: y        ! unknowns of the equations
        real(dp), dimension(:),     intent(out) :: dfdx     ! derivatives of f for x
        real(dp), dimension(:,:),   intent(out) :: dfdy     ! derivatives of f for y

        real(dp), dimension(:),     intent(in)  :: para      ! parameter

        real(dp) :: xtmp
        real(dp), dimension(size(y,1)) :: ytmp

        xtmp = x
        ytmp = y

        dfdx(1)     = 0.0_dp
        dfdx(2)     = 0.0_dp
        dfdx(3)     = 0.0_dp

        dfdy(1,1)   = -0.013_dp-para(1)*y(3)
        dfdy(1,2)   = 0.0_dp
        dfdy(1,3)   = -1000.0_dp*y(1)
        dfdy(2,1)   = 0.0_dp
        dfdy(2,2)   = -2500.0_dp*y(3)
        dfdy(2,3)   = -2500.0_dp*y(2)
        dfdy(3,1)   = -0.013_dp-1000.0_dp*y(3)
        dfdy(3,2)   = -2500.0_dp*y(3)
        dfdy(3,3)   = -1000.0_dp*y(1)-2500.0_dp*y(2)

    end subroutine jacobn2_testRB_para_dp

end module mo_eq_lv


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!     calculate the initial conditions for the gas
!         - the initial pressure
!         - identifying the chambers and gas pipe nodes
!           for gas networks
!         - the initial flow
!         - calculating the static temperature for gas networks
!
      subroutine initialnet(itg,ieg,ntg,ac,bc,lakon,v,
     &     ipkon,kon,nflow,ikboun,nboun,prop,ielprop,
     &     nactdog,ndirboun,nodeboun,xbounact,
     &     ielmat,ntmat_,shcon,nshcon,physcon,ipiv,nteq,
     &     rhcon,nrhcon,ipobody,ibody,xbodyact,co,nbody,network,
     &     iin_abs,vold,set,istep,iit,mi,ineighe,ilboun,channel,
     &     iaxial,nmpc,labmpc,ipompc,nodempc,coefmpc,ttime,time,
     &     iponoel,inoel)
!
      implicit none

      logical identity,gravity,gaspipe
!
      character*8 lakon(*)
      character*20 labmpc(*)
      character*81 set(*)
!
      integer mi(*),ieg(*),nflow,i,j,ntg,ielmat(mi(3),*),ntmat_,id,
     &     node1,node2,ider,iaxial,nmpc,nodempc(3,*),ipompc(*),
     &     nelem,index,nshcon(*),ipkon(*),kon(*),ikboun(*),nboun,idof,
     &     nodem,idirf(8),nactdog(0:3,*),imat,ielprop(*),id1,id2,
     &     nodef(8),ndirboun(*),nodeboun(*),itg(*),node,kflag,ipiv(*),
     &     nrhs,info,idof1,idof2,nteq,nrhcon(*),ipobody(2,*),ibody(3,*),
     &     nbody,numf,network,iin_abs,icase,index2,index1,nelem1,nelem2,
     &     node11,node21,node12,node22,istep,iit,ineighe(*),iponoel(*),
     &     ilboun(*),nelemup,k,node2up,idir,channel,inoel(2,*)
!
      real*8 ac(nteq,nteq), bc(nteq),prop(*),shcon(0:3,ntmat_,*),
     &     f,df(8),xflow,xbounact(*),v(0:mi(2),*),cp,r,tg1,
     &     tg2,gastemp,physcon(*),pressmin,dvi,rho,g(3),z1,z2,
     &     rhcon(0:1,ntmat_,*),co(3,*),xbodyact(7,*),kappa,
     &     a,Tt,Pt,Ts,pressmax,constant,vold(0:mi(2),*),href,
     &     coefmpc(*),ttime,time,xflow360
!
      kflag=1
      ider=0
      channel=0
      gaspipe=.false.
!
!     applying the boundary conditions
!
      do j=1,nboun
         v(ndirboun(j),nodeboun(j))=xbounact(j)
      enddo
!
!     check for channel elements
!
      do i=1,nflow
         nelem=ieg(i)
         if(lakon(nelem)(2:5).eq.'LICH') then
            channel=1
            return
         endif
      enddo
!
!     initializing ac and bc
!
      do i=1,nteq
         do j=1,nteq
            ac(i,j)=0.d0
         enddo
         bc(i)=0.d0
      enddo
!
!     for all but purely thermal networks:
!     determining the initial pressure
!     and identifying the chamber and gas pipe nodes
!
      if(network.gt.2) then
!
         call preinitialnet(ieg,lakon,v,ipkon,kon,nflow,prop,ielprop,
     &     ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,mi,iponoel,inoel,
     &     itg,ntg)
!
!     determining whether pressure initial conditions
!     are provided for all nodes
!
         pressmin=-1.d0
         pressmax=0.d0
         constant=1.55d0

         do i=1,ntg
            node=itg(i)
            if(v(2,node).lt.1.d-10) then
               v(2,node)=0.d0
            else
               if(pressmin.lt.0.d0) then
                  pressmin=v(2,node)
               elseif(v(2,node).lt.pressmin) then
                  pressmin=v(2,node)
               endif
!
               if(v(2,node).gt.pressmax)then
                  pressmax=v(2,node)
               endif
!
            endif
         enddo
!
         if(pressmin.lt.0.d0) then
            write(*,*)
     &           '*ERROR in initialnet: minimum initial pressure'
            write(*,*) '       is smaller than zero'
            call exit(201)
         endif
!
!     in nodes in which no initial pressure is given v(2,*)
!     is replaced by -n, where n is the number of elements the
!     node belongs to: allows to find boundary nodes of the
!     network
!
         do i=1,nflow
            nelem=ieg(i)
            index=ipkon(nelem)
            node1=kon(index+1)
            node2=kon(index+3)
            call nident(itg,node1,ntg,id1)
            call nident(itg,node2,ntg,id2)
!
            if (((lakon(nelem)(1:5).eq.'DGAPF').and.
     &           (iin_abs.eq.0))
     &           .or.((lakon(nelem)(1:3).eq.'DRE')
     &           .and.(lakon(nelem)(1:7).ne.'DREWAOR')
     &           .and.(iin_abs.eq.0))) then
!
!     In the case of a element of type GASPIPE or RESTRICTOR
!     (except TYPE= RESTRICTOR WALL ORIFICE)
!     the number of pipes connected to node 1 and 2
!     are computed and stored in ineighe(id1)
!     respectively ineighe(id2)
!
               gaspipe=.true.
               if(node1.ne.0) then
                  if (ineighe(id1).ge.0) then
!
                     if(node2.ne.0)then
                        ineighe(id1)=ineighe(id1)+1
                     endif
                  endif
               endif
               if(node2.ne.0) then
                  if (ineighe(id2).ge.0) then
                     if(node1.ne.0) then
                        ineighe(id2)=ineighe(id2)+1
                     endif
                  endif
               endif
            else
               if(iin_abs.eq.0) then
!
!     for all other elements (different from GASPIPE or
!     RESTRICTOR), including RESTRICTOR WALL ORIFICE
!     ineighe(idi)=-1
!     which means that they are connected to chambers
!     i.e. static and total values are equal
!
                  if (node1.ne.0) then
                     ineighe(id1)=-1
                  endif
                  if(node2.ne.0) then
                     ineighe(id2)=-1
                  endif
               endif
            endif
!
            if((node1.eq.0).or.(node2.eq.0)) cycle
            if(v(2,node1).lt.1.d-10) then
               v(2,node1)=v(2,node1)-1.d0
            endif
            if(v(2,node2).lt.1.d-10) then
               v(2,node2)=v(2,node2)-1.d0
            endif
         enddo
!
!     for each end node i: if ineighe(i)<0: chamber
!     else: ineighe(i)=number of pipe connections
!
!     assigning values to the boundary nodes of the network
!     (i.e. nodes belonging to only one element)
!
         do i=1,ntg
            node=itg(i)
!     boundary condition nodes
            if(nactdog(2,node).eq.0) then
               v(2,node)=dtan(v(2,node)
     &              *constant/pressmax)
               cycle
            endif
!     initial condition nodes
            if((nactdog(2,node).ne.0)
     &           .and.(v(2,node).gt.0d0)) then
               v(2,node) = dtan(v(2,node)
     &              *constant/pressmax)
               cycle
            endif
!     nodes neither defined as *BOUNDARY nor as *INITIAL CONDITIONS
            if(abs(v(2,node)+1.d0).lt.1.d-10) then
               v(2,node)=0.95d0*pressmin
               pressmin=0.95d0*pressmin
               v(2,node)=dtan(v(2,node)
     &              *constant/pressmax)
            endif
         enddo
!
!     mass flow and geometry dofs: 1 on the diagonal
!
         do i=1,ntg
            node=itg(i)
            do j=1,1
               if(nactdog(j,node).ne.0) then
                  idof=nactdog(j,node)
                  ac(idof,idof)=1.d0
               endif
            enddo
            do j=3,3
               if(nactdog(j,node).ne.0) then
                  idof=nactdog(j,node)
                  ac(idof,idof)=1.d0
               endif
            enddo
         enddo
!
!     pressure dofs
!
         do i=1,nflow
            nelem=ieg(i)
            index=ipkon(nelem)
            node1=kon(index+1)
            node2=kon(index+3)
            if((node1.eq.0).or.(node2.eq.0)) cycle
            idof1=nactdog(2,node1)
            idof2=nactdog(2,node2)
            if(idof1.ne.0) then
               if(v(2,node1).gt.0.d0) then
!     initial pressure given in node1
                  ac(idof1,idof1)=1.d0
                  bc(idof1)=v(2,node1)
               else
                  ac(idof1,idof1)=ac(idof1,idof1)+1.d0
                  if(v(2,node2).gt.0.d0) then
!     initial pressure given in node2
                     bc(idof1)=bc(idof1)+v(2,node2)
                  else
                     ac(idof1,idof2)=ac(idof1,idof2)-1.d0
                  endif
               endif
            endif
            if(idof2.ne.0) then
               if(v(2,node2).gt.0.d0) then
!     initial pressure given in node2
                  ac(idof2,idof2)=1.d0
                  bc(idof2)=v(2,node2)
               else
                  ac(idof2,idof2)=ac(idof2,idof2)+1.d0
                  if(v(2,node1).gt.0.d0) then
!     initial pressure given in node1
                     bc(idof2)=bc(idof2)+v(2,node1)
                  else
                     ac(idof2,idof1)=ac(idof2,idof1)-1.d0
                  endif
               endif
            endif
         enddo
!
!     checking for nodes not belonging to network elements
!
         do i=1,ntg
            node=itg(i)
            if(nactdog(2,node).ne.0) then
               idof=nactdog(2,node)
               if(ac(idof,idof).eq.0.d0) then
                  ac(idof,idof)=1.d0
                  bc(idof)=v(2,node)
               endif
            endif
         enddo
      endif
!
!     temperature conditions for all but purely hydrodynamic
!     networks
!
      if(network.ne.4) then
!
!     temperature dofs
!
         do i=1,nflow
            nelem=ieg(i)
            index=ipkon(nelem)
            node1=kon(index+1)
            nodem=kon(index+2)
            node2=kon(index+3)
            if((node1.eq.0).or.(node2.eq.0)) cycle
            idof1=nactdog(0,node1)
            idof2=nactdog(0,node2)
            if(idof1.ne.0) then
               if(v(0,node1).gt.0.d0) then
!     initial temperature given in node1
                  ac(idof1,idof1)=1.d0
                  bc(idof1)=v(0,node1)
!
!                 temperature taken into account only for incoming
!                 flux (from the viewpoint of node1).
!
               elseif(v(1,nodem).lt.0.d0) then
                  ac(idof1,idof1)=ac(idof1,idof1)+1.d0
                  if(v(0,node2).gt.0.d0) then
!     initial temperature given in node2
                     bc(idof1)=bc(idof1)+v(0,node2)
                  else
                     ac(idof1,idof2)=ac(idof1,idof2)-1.d0
                  endif
               endif
            endif
            if(idof2.ne.0) then
               if(v(0,node2).gt.0.d0) then
!     initial temperature given in node2
                  ac(idof2,idof2)=1.d0
                  bc(idof2)=v(0,node2)
!
!                 temperature taken into account only for incoming
!                 flux (from the viewpoint of node2). Default (if
!                 no flux is predefined) is positive flux.
!
               elseif(v(1,nodem).ge.0.d0) then
                  ac(idof2,idof2)=ac(idof2,idof2)+1.d0
                  if(v(0,node1).gt.0.d0) then
!     initial temperature given in node1
                     bc(idof2)=bc(idof2)+v(0,node1)
                  else
                     ac(idof2,idof1)=ac(idof2,idof1)-1.d0
                  endif
               endif
            endif
         enddo
!
!     checking for nodes not belonging to network elements
!
         do i=1,ntg
            node=itg(i)
            if(nactdog(0,node).ne.0) then
               idof=nactdog(0,node)
               if(ac(idof,idof).eq.0.d0) then
                  ac(idof,idof)=1.d0
                  bc(idof)=v(0,node)
               endif
            endif
         enddo
      endif
!
!     additional multiple point constraints
!
      do j=1,nteq
         if(dabs(ac(j,j)).lt.1.d-20) ac(j,j)=1.d0
      enddo
!
!     solving the system
!
      if(nteq.gt.0) then
         nrhs=1
         call dgesv(nteq,nrhs,ac,nteq,ipiv,bc,nteq,info)
         if(info.ne.0) then
            write(*,*) '*ERROR in initialnet: singular matrix'
            call exit(201)
         endif
      endif
!
!     storing the initial pressure in v
!     (not for purely thermal networks)
!
      if(network.gt.2) then
         do i=1,ntg
            node=itg(i)
            if(nactdog(2,node).eq.0) then
               v(2,node)=pressmax
     &              *datan(v(2,node))/constant
               cycle
            endif
            v(2,node)=pressmax
     &           *datan(bc(nactdog(2,node)))/constant
         enddo
      endif
!
!     storing the initial temperature in v
!     (not for purely hydrodynamic networks)
!
      if(network.ne.4) then
         do i=1,ntg
            node=itg(i)
            if(nactdog(0,node).ne.0) v(0,node)=bc(nactdog(0,node))
         enddo
      endif
!
!     for ELEMENT TYPE BRANCH
!     initial pressures in branch 1 and 2 are equal and set to the
!     maximum of the two pressures
!
      do i=1,nflow
         nelem=ieg(i)
         if(lakon(i)(2:5).ne.'REBR') cycle
         index=ielprop(nelem)
!
         nelem1=nint(prop(index+2))
         index1=ipkon(nelem1)
         node11=kon(index1+1)
         node21=kon(index1+3)
!
         nelem2=nint(prop(index+3))
         index2=ipkon(nelem2)
         node12=kon(index2+1)
         node22=kon(index2+3)
!
         if(node11.eq.node12) then
            node1=node21
            node2=node22
         elseif(node11.eq.node22) then
            node1=node21
            node2=node12
         elseif(node21.eq.node12) then
            node1=node11
            node2=node22
         elseif(node21.eq.node22) then
            node1=node11
            node2=node12
         endif
!
         if(lakon(i)(6:6).eq.'S') then
!
!           maximum
!
            if(v(2,node1).ge.v(2,node2)) then
               v(2,node2)=v(2,node1)
            else
               v(2,node1)=v(2,node2)
            endif
         elseif(lakon(i)(6:6).eq.'J') then
!
!           minimum
!
            if(v(2,node1).ge.v(2,node2)) then
               v(2,node1)=v(2,node2)
            else
               v(2,node2)=v(2,node1)
            endif
         endif
      enddo
!
!     identifying the chamber nodes for purely thermal
!     gas networks (needed to determine the static
!     temperature which is used for the material properties)
!
      if(network.le.2) then
         do i=1,nflow
            nelem=ieg(i)
            if((lakon(nelem)(2:3).eq.'LP').or.
     &           (lakon(nelem)(2:3).eq.'LI')) cycle
            index=ipkon(nelem)
            node1=kon(index+1)
            node2=kon(index+3)
            if(node1.ne.0) then
               call nident(itg,node1,ntg,id1)
               ineighe(id1)=-1
            endif
            if(node2.ne.0) then
               call nident(itg,node2,ntg,id2)
               ineighe(id2)=-1
            endif
         enddo
      endif
!
!     initialisation of bc
!
      do i=1,nteq
         bc(i)=0.d0
      enddo
!
!     determining the initial mass flow in those nodes for which no
!     flux boundary conditions are defined
!
      if(network.gt.2)then
         do i=1,nflow
            nelem=ieg(i)
!
            index=ipkon(nelem)
            node1=kon(index+1)
            node2=kon(index+3)
            if((node1.eq.0).or.(node2.eq.0)) cycle
            nodem=kon(index+2)
!
!     cycle if both the geometry and the mass flow are known
!
            if((nactdog(1,nodem).eq.0).and.(nactdog(3,nodem).eq.0))
     &           cycle
!
!     for liquids: determine the gravity vector
!
            if(lakon(nelem)(2:3).eq.'LI') then
               gravity=.false.
               do j=1,3
                  g(j)=0.d0
               enddo
               if(nbody.gt.0) then
                  index=nelem
                  do
                     j=ipobody(1,index)
                     if(j.eq.0) exit
                     if(ibody(1,j).eq.2) then
                        g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j)
                        g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j)
                        g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j)
                        z1=-g(1)*co(1,node1)-g(2)*co(2,node1)
     &                       -g(3)*co(3,node1)
                        z2=-g(1)*co(1,node2)-g(2)*co(2,node2)
     &                       -g(3)*co(3,node2)
                        gravity=.true.
                     endif
                     index=ipobody(2,index)
                     if(index.eq.0) exit
                  enddo
               endif
               if(.not.gravity) then
                  write(*,*)
     &                 '*ERROR in initialnet: no gravity vector'
                  write(*,*)
     &                 '       was defined for liquid element',nelem
                  call exit(201)
               endif
            endif
!
            tg1=v(0,node1)
            tg2=v(0,node2)
!
!     For liquid pipes elements the upstream temperature
!     is used to determine the material properties
!     For all other elements the averaged value of the
!     temperature between inlet and outlet is used
!
            if((lakon(nelem)(2:3).ne.'LP').and.
     &           (lakon(nelem)(2:3).ne.'LI')) then
               gastemp=(tg1+tg2)/2.d0
            else
               xflow=v(1,nodem)
               if(xflow.gt.0) then
                  gastemp=tg1
               else
                  gastemp=tg2
               endif
            endif
!
            imat=ielmat(1,nelem)
            call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,cp,
     &           r,dvi,rhcon,nrhcon,rho)
!
!     If inlet and outlet pressure are equal this leads to a massflow rate
!     equal to 0 and in turn possibly to a singular matrix configuration
!
!     gas
            if(((lakon(nelem)(2:3).ne.'LI').and.
     &           (v(2,node1).eq.v(2,node2))).or.
!     liquid
     &           ((lakon(nelem)(2:3).eq.'LI').and.
     &           (v(2,node1)+z1.eq.v(2,node2)+z2))) then

!
!     if neither inlet nor outlet pressure are active D.O.F: error-> stop
!
               if((nactdog(2,node1).eq.0)
     &           .and.(nactdog(2,node2).eq.0)) then
                  WRITE(*,*) '**************************************'
                  write(*,*) '*ERROR:in subroutine initialnet.f'
                  write(*,*) '       in element', nelem
                  write(*,*) '       Inlet and outlet pressures are '
                  write(*,*) '       boundary conditions '
                  write(*,*) '       node1',node1,' pressure',
     &                 v(2,node1)
                  write(*,*) '       node2',node2,' pressure',
     &                 v(2,node2)
                  call exit(201)
!
!     if inlet pressure is an active degree of freedom
!
               else if((nactdog(2,node1).ne.0)
     &                 .and.(nactdog(2,node2).eq.0))then
                  WRITE(*,*) '**************************************'
                  write(*,*) '*WARNING:in subroutine initialnet.f'
                  write(*,*) '       in element', nelem
                  write(*,*)
     &                 '       Inlet pressure initial condition '
                  write(*,*) '       is changed '
                  write(*,*) '       node1',node1,
     &                 ' given initial pressure',v(2,node1)
                  v(2,node1)=1.1*v(2,node1)
                  write(*,*) '       node1',node1,
     &                 ' new initial pressure',v(2,node1)
                  write(*,*) '       node2',node2,' pressure',
     &                 v(2,node2)
!
!     if outlet pressure is an active D.O.F.
!
               else if((nactdog(2,node1).eq.0)
     &                 .and.(nactdog(2,node2).ne.0))then
                  WRITE(*,*) '**************************************'
                  write(*,*) '*WARNING:in subroutine initialnet.f'
                  write(*,*) '       in element', nelem
                  write(*,*)
     &                 '       Outlet pressure initial condition '
                  write(*,*) '       is changed '
                  write(*,*) '       node1',node1,' pressure'
     &                 ,v(2,node1)
                  write(*,*) '       node2',node2,
     &                 'given intial pressure',
     &                 v(2,node2)
                  v(2,node2)=0.9*v(2,node2)
                  write(*,*) '       node2',node2,
     &                 ' new initial pressure',v(2,node2)
!
!     if both inlet and outlet pressures are active D.O.F.
!
               else if((nactdog(2,node1).ne.0)
     &                 .and.(nactdog(2,node2).ne.0))then
                  WRITE(*,*) '**************************************'
                  write(*,*) '*WARNING:in subroutine initialnet.f'
                  write(*,*) '       in element', nelem
                  write(*,*) '       Inlet and outlet pressure '
                  write(*,*) '       initial condition are changed '
                  write(*,*) '       node1',node1,
     &                 ' given initial pressure',v(2,node1)
                  v(2,node1)=1.05*v(2,node2)
                  write(*,*) '       node1',node1,
     &                 ' new intial pressure',v(2,node1)
                  write(*,*) '       node2',node2,
     &                 ' given initial pressure',v(2,node2)
                  v(2,node2)=0.95*v(2,node2)
                  write(*,*) '       node2',node2,
     &                 ' new intial pressure',v(2,node2)
               endif
            endif
!
!           calculating flux if the flux is an unknown AND there was
!           no initial flux defined by the user
!
            if((nactdog(1,nodem).ne.0).and.(v(1,nodem).eq.0.d0)) then
               call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
     &           nactdog,identity,ielprop,prop,kflag,v,xflow,f,
     &           nodef,idirf,df,cp,r,rho,physcon,g,co,dvi,numf,
     &           vold,set,shcon,nshcon,rhcon,nrhcon,ntmat_,mi,ider,
     &           ttime,time,iaxial)
               v(1,nodem)=xflow
            endif
!
            if((lakon(nelem)(2:4).ne.'LIP').and.
     &         (lakon(nelem)(2:3).ne.'VO').and.
     &         (lakon(nelem)(2:4).ne.'ATR').and.
     &         (lakon(nelem)(2:4).ne.'RTA')) then
               if(v(1,nodem).eq.0d0) then
                  WRITE(*,*) '**************************************'
                  write(*,*) '*ERROR:in subroutine initialnet.f'
                  write(*,*) '       in element', nelem,
     &                 lakon(nelem)(1:6)
                  write(*,*) '       mass flow rate value = 0 !'
                  write(*,*) '       node1',node1,' pressure',
     &                 v(2,node1)
                  write(*,*) '       node2',node2,' pressure',
     &                 v(2,node2)
                  call exit(201)
               endif
               if (v(1,nodem)*(v(2,node1)-v(2,node2)).lt.0) then
                  WRITE(*,*) '**************************************'
                  write(*,*) '*WARNING: in subroutine initialnet.f'
                  write(*,*) '        in element', nelem
                  write(*,*)
     &             '        initial mass flow rate value does not'
                  write(*,*) '        correspond to pressure gradient'
                  write(*,*) '        node1',node1,'pressure',
     &                 v(2,node1)
                  write(*,*) '        node2',node2,'pressure',
     &                 v(2,node2)
                  write(*,*) '        nodem',nodem,'mass flow',
     &                 v(1,nodem)
                  write(*,*) '        the sign of the initial mass flow'
                  write(*,*) '        rate will be changed'
c                  v(1,nodem)=-v(1,nodem)
               endif
c               if (v(1,nodem).lt.0) then
c                  WRITE(*,*) '**************************************'
c                  write(*,*) '*WARNING: in subroutine initialnet.f'
c                  write(*,*) '        in element', nelem
c                  write(*,*) '        mass flow rate value .le. 0 !'
c                  write(*,*) '        node1',node1,'pressure',
c     &                 v(2,node1)
c                  write(*,*) '        node2',node2,'pressure',
c     &                 v(2,node2)
c                  write(*,*) '        check element definition'
c               endif
            endif
         enddo
      endif
!
      call postinitialnet(ieg,lakon,v,ipkon,kon,nflow,prop,ielprop,
     &     ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,mi,iponoel,inoel,
     &     itg,ntg)
!
!     calculating the static temperature for nodes belonging to gas pipes
!     and restrictors (except RESTRICTOR WALL ORIFICE)
!
      if (gaspipe.and.(iin_abs.eq.0)) then
!
!     ineighe(i) is set to -1 (= chamber node) if more than 2 pipes
!     are connected to the node
!
         do i=1,ntg
            if(ineighe(i).gt.2) then
               ineighe(i)=-1
               write(*,*) '*WARNING :in subroutine initialnet.f'
               write(*,*) '          more than 2 elements GASPIPE'
               write(*,*) '          or RESTRICTOR are connected '
               write(*,*) '          to node',itg(i),'. The common'
               write(*,*)
     &              '          node is converted into a chamber.'
               write(*,*) '          Total and static parameters are'
               write(*,*) '          equal'
            endif
         enddo
!
         do i=1,nflow
            nelem=ieg(i)
            index=ipkon(nelem)
            node1=kon(index+1)
            node2=kon(index+3)
            call nident(itg,node1,ntg,id1)
            call nident(itg,node2,ntg,id2)
!
!           for each end node i:
!             if ineighe(i)=-1: chamber
!                          = 0: middle node
!                          > 0: number of pipe connections (max. 2 allowed)
!
!           the number of pipe connections for ineighe(i)>0 is
!           replaced by the number of an element the node belongs to:
!
            if(node1.gt.0) then
               if((ineighe(id1).eq.1).or.
     &              (ineighe(id1).eq.2)) then
                  if(node2.ne.0)then
                     ineighe(id1)=nelem
                  endif
               endif
            endif
!
            if(node2.gt.0) then
               if((ineighe(id2).eq.1).or.
     &              (ineighe(id2).eq.2)) then
                  if (node1.ne.0) then
                     ineighe(id2)=nelem
                  endif
               endif
            endif
         enddo
!
!     The static temperature is calculated and stored in v(3,node)
!     total temperatures are supposed equal (adiabatic pipe)
!
         do i=1,ntg
            node=itg(i)
            if(ineighe(i).gt.0) then
!
               nelem=ineighe(i)
               index=ielprop(nelem)
               nodem=kon(ipkon(nelem)+2)
!
               imat=ielmat(1,nelem)
               call materialdata_tg(imat,ntmat_,v(0,node),
     &              shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho)
               kappa=cp/(cp-R)
               xflow=v(1,nodem)
               Tt=v(0,node)
               Pt=v(2,node)
!
               if(lakon(nelem)(2:5).eq.'GAPF') then
                  A=prop(index+1)
                  if(lakon(nelem)(2:6).eq.'GAPFA') then
                     icase=0
                  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
                     icase=1
                  endif
!
               elseif(lakon(nelem)(2:3).eq.'RE') then
                  index2=ipkon(nelem)
                  node1=kon(index2+1)
                  node2=kon(index2+3)
                  if(lakon(nelem)(4:5).eq.'EX') then
                     if(lakon(nint(prop(index+4)))(2:6).eq.'GAPFA')
     &                    then
                        icase=0
                     elseif
     &                   (lakon(nint(prop(index+4)))(2:6).eq.'GAPFI')
     &                       then
                        icase=1
                     endif
                  else
                     icase=0
                  endif
!
                  if(lakon(nelem)(4:5).eq.'BE') then
                     a=prop(index+1)
!
                  elseif(lakon(nelem)(4:5).eq.'BR') then
                     if(lakon(nelem)(4:6).eq.'BRJ') then
                        if(nelem.eq.nint(prop(index+2)))then
                           A=prop(index+5)
                        elseif(nelem.eq.nint(prop(index+3))) then
                           A=prop(index+6)
                        endif
                     elseif(lakon(nelem)(4:6).eq.'BRS') then
                        if(nelem.eq.nint(prop(index+2)))then
                           A=prop(index+5)
                        elseif(nelem.eq.nint(prop(index+3))) then
                           A=prop(index+6)
                        endif
                     endif
!
                  else
                     if(node.eq.node1) then
                        a=prop(index+1)
                     elseif(node.eq.node2) then
                        a=prop(index+2)
                     endif
                  endif
               endif
!
               if(v(3,node).eq.0.d0) then
                  xflow360=xflow*iaxial
                  call ts_calc(xflow360,Tt,Pt,kappa,r,a,Ts,icase)
                  v(3,node)=Ts
               endif
!
!     if the element is not of gaspipe or branch type,
!     total and static temperatures are equal for all endnodes
!
            endif
         enddo
      endif
!
!     for chambers the static temperature equals the total
!     temperature
!
      do i=1,ntg
         if(ineighe(i).eq.-1) v(3,itg(i))=v(0,itg(i))
c         write(*,*) 'initialnet ',i,v(0,itg(i)),
c     &      v(1,itg(i)),v(2,itg(i)),ineighe(i)
      enddo
!
c      write(*,*) 'initialnet '
c      do i=1,ntg
c         write(*,'(i10,3(1x,e11.4))') itg(i),(v(j,itg(i)),j=0,2)
c      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
C------------------------------------------------------------------
C Get navigation data for reference image.
C
      SUBROUTINE GET_REF_NAV(PROJECT,TARGET,ITYPE,sbuf2,ibuf2,t2)
      IMPLICIT NONE
      CHARACTER*5 PROJECT
      CHARACTER*12 TARGET
      REAL*8 SBUF2(100)
      INTEGER*4 ITYPE,IBUF2(200)
      REAL*8 T2

      COMMON/C2/OM2(3,3),RS2(3),RX2(3),RX2MAG
      COMMON/C2/FL2,OAL2,OAS2,SCALE2,ZSCALE2
      COMMON/C2/ICAM2,ITYPE2,CONV2(3),PROJECT2
      REAL*8 OM2,RS2,RX2,RX2MAG,FL2,OAL2,OAS2,SCALE2,ZSCALE2
      INTEGER*4 ICAM2,ITYPE2,CONV2
      CHARACTER*5 PROJECT2

      INTEGER*4 LAB(80),DATA(40),DEF,RUNIT,NUM,IND,IFDSC,I,ISOURCE
      REAL*4 XFL,XOAL,XOAS,XSCALE
      REAL*8 MP

      CHARACTER*128 REFFILE

C     ....Get info for the reference image (if given)
      CALL XVPARM('REF',reffile,num,def,1)
      IF (NUM.NE.1) GOTO 10
      CALL XVUNIT(runit,'N/A',1,ind,'U_NAME',REFFILE,' ')
      CALL XVOPEN(RUNIT,ind, 'OPEN_ACT', 'SA', 'IO_ACT','SA',' ')
      CALL GETPROJ(RUNIT,project2,icam2,ifdsc,ind)
      IF (IND.NE.0) CALL MABEND('***Invalid ref image project ID')
      CALL XVMESSAGE('Reference image spacecraft is '//PROJECT2,' ')
      CALL GETLABCON(RUNIT,project2,lab,ind)
      IF (IND.GT.1) CALL MABEND('***Err reading ref image label')
      ICAM2 = LAB(6)

C     ...Determine if input image is geometrically corrected
      CALL MP_INIT(mp,ind)
      CALL MP_LABEL_READ(MP,RUNIT,ind)
      IF (IND.GE.0) THEN     !Map labels found
         CALL MP_MPO2BUF(MP,data,ind)
         ITYPE = DATA(39)
      ELSE
         CALL SEARC_DISTOR(RUNIT,ind)
         ITYPE2 = 7             !Geometrically uncorrected (image-space)
         IF (IND.EQ.1) ITYPE2=8	!Geometrically corrected (object_space)
      ENDIF
      IF (ITYPE2.EQ.7 .AND. (PROJECT2.EQ.'VIKOR'.OR.
     +          PROJECT2.EQ.'VGR-1'.OR.PROJECT2.EQ.'VGR-2'))
     +  CALL MABEND('Reference image must be geometricaly corrected')
      GOTO 30
C
C     ....Here if reference frame is not specified
   10 CALL XVMESSAGE('No reference image provided',' ')
      DO I=1,80
         LAB(I) = -999		!Initialize dummy label buffer
      ENDDO
      CALL MVCL(TARGET,lab(25),12)
      CALL XVPARM('RMISSION',project2,num,def,1)
      IF (DEF.EQ.1) PROJECT2=PROJECT
      CALL XVPARM('RSCET',lab(8),num,def,6)
      IF (DEF.EQ.1) CALL MABEND('***RSCET parameter required')
      CALL XVPARM('RCAM',icam2,num,def,1)	!reference camera S/N
      IF (DEF.EQ.1) CALL MABEND('***RCAM parameter required')
      LAB(6) = ICAM2
      ITYPE2 = ITYPE

C     ....Get camera focal length, optical axis line,samp, picture scale
   30 CALL GETCAMCON(PROJECT2,ICAM2,xfl,xoal,xoas,xscale,ind)
      IF (IND.NE.0) CALL MABEND('***Invalid camera S/N for ref image')
      FL2 = XFL
      OAL2 = XOAL
      OAS2 = XOAS
      SCALE2 = XSCALE
      ZSCALE2 = FL2*SCALE2
      IF (ITYPE2.EQ.7) THEN
         CALL MVCL(PROJECT2,conv2,5)
         CONV2(3) = ICAM2
      ENDIF

C     ....Get geometry data for reference image
      CALL GETSPICE4(PROJECT2,.TRUE.,LAB,sbuf2,ind)
      IF (IND.NE.1) CALL MABEND('***Err getting SPICE for ref image')
      CALL CMSOURCE(SBUF2,isource)	!Print C matrix source
      CALL MVE(8,9,SBUF2(59),om2,1,1)		!OM matrix
      CALL MVE(8,3,SBUF2(22),rs2,1,1)		!RS vector
      CALL GET_TIME(SBUF2,t2)
      RETURN
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
The (Western Pacific) Depot is a gay bar located on the corner of 20th and K Street, directly across from Faces. They have 30+ TV screens for sports and music videos, an enclosed patio in the back (where you can smoke), and a pool table. They have up to 4 bartenders on Friday and Saturday nights. Everyone usually floods the front bar, so try heading to the back to get your drinks. On Wednesday nights they have a pretty popular pop culture trivia quiz.

Drinks are pulled off the floor around 1:45 AM. (Its illegal to have booze on the floor past 2 AM in Sacramento.) But they stay open until 4 AM on Fridays and Saturdays to serve nonalcoholic beverages.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program ufsLandDriver

  use ufsLandNoahDriverModule
  use ufsLandNoahMPDriverModule
  use ufsLandNoahType, only      : noah_type
  use ufsLandNoahMPType, only    : noahmp_type
  use NamelistRead
  use ufsLandStaticModule, only  : static_type
  use ufsLandForcingModule, only : forcing_type

  implicit none

  type (noah_type)     :: noah
  type (noahmp_type)   :: noahmp
  type (namelist_type) :: namelist
  type (static_type)   :: static
  type (forcing_type)  :: forcing

  integer, parameter :: NOAH_LAND_SURFACE_MODEL = 1
  integer, parameter :: NOAHMP_LAND_SURFACE_MODEL = 2

  call namelist%ReadNamelist()

  namelist%begsub = namelist%begloc
  namelist%endsub = namelist%endloc
  namelist%lensub = namelist%endloc - namelist%begloc + 1

  land_model : select case(namelist%land_model)

    case(NOAH_LAND_SURFACE_MODEL)

      call ufsLandNoahDriverInit(namelist, static, forcing, noah)

      call ufsLandNoahDriverRun(namelist, static, forcing, noah)

      call ufsLandNoahDriverFinalize()

    case(NOAHMP_LAND_SURFACE_MODEL)

      call ufsLandNoahMPDriverInit(namelist, static, forcing, noahmp)

      call ufsLandNoahMPDriverRun(namelist, static, forcing, noahmp)

      call ufsLandNoahMPDriverFinalize()

    case default

      stop "no valid land_model set in namelist"

  end select land_model

end program

max_stars_repo_name: JCSDA-internal/UCLDAS

max_issues_repo_name: JCSDA-internal/UCLDAS

max_forks_repo_name: JCSDA-internal/ucldas


text:
http://mortgageindavis.com/

Marc Thompson  Mortgage Planner.  Powered by UmaxMortgage.

Check out Marcs Facebook page at SaveMoneyEnjoyLife.  Lots of interesting trivia, pics.

Marc Thompson, mthompson@umaxmortgage.com; (530) 7538800.   437 F Street, Davis, CA MonFri 8:305:30.

Need a new mortgage, refinance, reverse mortgage?  Marc can help with all of those.  Great personal service.  Give him a call.


text:
program my_fib
        use mytest
        integer :: N = 8
        !$omp declare target(fib)
        !$omp target map(tofrom: N)
        call fib(N)
        !$omp end target
        write(*, *) "N=", N
end program

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
!!# MODULE <<FUN_Reorder>>
MODULE FUN_Reorder

!!## PURPOSE
!! Performs the a reordering of a list according to the argument "order",
!! which gives the new location of each element.  Useful for reordering
!! something else according to the same order of a Sort operation.  Ignores
!! values of order <=0.

!!## USAGE
!
!      new_list = Reorder( old_list , order [, side ])
!
!! where <old_list> is a 1D array of any intrinsic type,
!! <order> is a 1D array of  the same size as <old_list>.
!!
!! On return, < new_list(order(i)) = old_list(i) >,
!! forall <order(i)>0>, in the default case or if <side="L">.
!! In the case of <side="R">, right-side reordering is performed
!! which leads to a different output,  < new_list(i) = old_list(order(i)) >,
!! forall <order(i)>0>.

!!## NOTES
!! The return value of the string reordering <Reorder_S> is the
!! variable length string <varying_string>.  I wish it could be the
!! a normal string with length set by <LEN> but the compiler
!! will not allow it.

!!## MODULES
USE KND_IntrinsicTypes,ONLY: KIND_Rsp,KIND_Rdp,& !!((01-A-KND_IntrinsicTypes.f90))
              KIND_Csp,KIND_Cdp,&
              KIND_I1,KIND_I2,KIND_I4,KIND_I8,&
              KIND_L1,KIND_L2,KIND_L4,&
              KIND_S
USE ISO_varying_string                           !!((03-A-ISO_varying_string.f90))
USE FUN_Default                                  !!((04-A-FUN_Default.f90))

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## DEFAULT ACCESS
PRIVATE

!!## INTERFACES
INTERFACE Reorder
 MODULE PROCEDURE Reorder_Rsp
 MODULE PROCEDURE Reorder_Rdp
 MODULE PROCEDURE Reorder_I1
 MODULE PROCEDURE Reorder_I2
 MODULE PROCEDURE Reorder_I4
 MODULE PROCEDURE Reorder_I8
 MODULE PROCEDURE Reorder_L1
 MODULE PROCEDURE Reorder_L2
 MODULE PROCEDURE Reorder_L4
 MODULE PROCEDURE Reorder_S
 MODULE PROCEDURE Reorder_VS
END INTERFACE

!!## PUBLIC ACCESS LIST
PUBLIC :: Reorder

!!## MODULE PROCEDURES
CONTAINS

!!### PURE FUNCTION  <<Reorder_Rsp>>
PURE FUNCTION Reorder_Rsp( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rsp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_R.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_Rdp>>
PURE FUNCTION Reorder_Rdp( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rdp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_R.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_Csp>>
PURE FUNCTION Reorder_Csp( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_C=>KIND_Csp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_C.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_Cdp>>
PURE FUNCTION Reorder_Cdp( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_C=>KIND_Cdp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_C.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_I1>>
PURE FUNCTION Reorder_I1( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I1     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_I.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_I2>>
PURE FUNCTION Reorder_I2( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I2     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_I.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_I4>>
PURE FUNCTION Reorder_I4( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I4     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_I.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_I8>>
PURE FUNCTION Reorder_I8( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I8     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_I.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_L1>>
PURE FUNCTION Reorder_L1( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L1     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_L.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_L2>>
PURE FUNCTION Reorder_L2( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L2     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_L.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_L4>>
PURE FUNCTION Reorder_L4( old_list , order , side )  RESULT(new_list)
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L4     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-A-FUN_Reorder_L.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_S>>
PURE FUNCTION Reorder_S( old_list , order , side )  RESULT(new_list)
INCLUDE "05-A-FUN_Reorder_S.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

!!### PURE FUNCTION  <<Reorder_VS>>
PURE FUNCTION Reorder_VS( old_list , order , side )  RESULT(new_list)
INCLUDE "05-A-FUN_Reorder_VS.f90.hdr"
!!--begin--
INCLUDE "05-A-FUN_Reorder.f90.bdy"
!!--end--
END FUNCTION

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack