text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!--------------------------------------------------------------------
!
!     Some useful debugging subroutines
!
!--------------------------------------------------------------------

!     Writes solution data and residue to file for debugging. Can also
!     be used with multiple processes
      SUBROUTINE DEBUGR()
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      INTEGER(KIND=IKIND) a, i, Ac, fid, iM, nNo, gnNo, s, e, ldof
      CHARACTER(LEN=stdL) fName

      REAL(KIND=RKIND), ALLOCATABLE :: lX(:,:), lA(:,:), lY(:,:),
     2   lD(:,:), lR(:,:), gX(:,:), gA(:,:), gY(:,:), gD(:,:), gR(:,:)

      fid = 1889
      WRITE(fName,'(A)') TRIM(appPath)//"dbgR_"//TRIM(eq(cEq)%sym)//
     2   "_"//STR(cTS)//"_"//STR(eq(cEq)%itr)
      IF (cm%mas()) OPEN(fid,FILE=TRIM(fName))

      s = eq(cEq)%s
      e = eq(cEq)%e
      ldof = e-s+1
      DO iM=1, nMsh
         nNo  = msh(iM)%nNo
         gnNo = msh(iM)%gnNo
         ALLOCATE(lX(nsd,nNo), lA(lDof,nNo), lY(lDof,nNo), lD(lDof,nNo),
     2      lR(dof,nNo))
         IF (cm%mas()) THEN
            ALLOCATE(gX(nsd,gnNo), gA(lDof,gnNo), gY(lDof,gnNo),
     2       gD(lDof,tnNo), gR(dof,gnNo))
         ELSE
            ALLOCATE(gX(0,0), gA(0,0), gY(0,0), gD(0,0), gR(0,0))
         END IF
         lX = 0._RKIND
         lA = 0._RKIND
         lY = 0._RKIND
         lD = 0._RKIND
         lR = 0._RKIND
         DO a=1, msh(iM)%nNo
            Ac = msh(iM)%gN(a)
            lX(:,a) = x(:,Ac)
            lA(:,a) = An(s:e,Ac)
            lY(:,a) = Yn(s:e,Ac)
            lD(:,a) = Dn(s:e,Ac)
            lR(:,a) = R(:,Ac)
         END DO

         gX = GLOBAL(msh(iM), lX)
         gA = GLOBAL(msh(iM), lA)
         gY = GLOBAL(msh(iM), lY)
         gD = GLOBAL(msh(iM), lD)
         gR = GLOBAL(msh(iM), lR)

         IF (cm%mas()) THEN
            WRITE(fid,'(A)') "Mesh: <"//TRIM(msh(iM)%name)//">"
            DO a=1, gnNo
               WRITE(fid,'(A)',ADVANCE='NO') STR(a)
               DO i=1, nsd
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gX(i,a))
               END DO
               DO i=1, dof
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gR(i,a))
               END DO
               DO i=1, lDof
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gA(i,a))
               END DO
               DO i=1, lDof
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gY(i,a))
               END DO
               DO i=1, lDof
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gD(i,a))
               END DO
               WRITE(fid,'(A)')
            END DO
         END IF
         DEALLOCATE(lX, lA, lY, lD, lR, gX, gA, gY, gD, gR)
      END DO

      IF (cm%mas()) CLOSE(fid)

      iM = cm%reduce(iM)

      RETURN
      END SUBROUTINE DEBUGR
!--------------------------------------------------------------------
      SUBROUTINE PDEBUGVALR()
      USE COMMOD
      IMPLICIT NONE

      INTEGER(KIND=IKIND) a, b, i, j, k, fid
      CHARACTER(LEN=stdL) fName

      fid = 1256 + cm%tF()
      i   = eq(cEq)%itr
      WRITE(fName,'(A)') TRIM(appPath)//"Val_R_"//STR(cTS)//"_"//
     2   STR(i)//"_"//STR(cm%tF())
      OPEN(fid,FILE=TRIM(fName))
      DO a=1, tnNo
         WRITE(fid,'(A)') REPEAT('=',32)
         WRITE(fid,'(A)',ADVANCE='NO') "Row: "//STR(a)//" grow: "//
     2      STR(ltg(a))//" x: "

         DO i=1, nsd
            WRITE(fid,'(A)',ADVANCE='NO') " "//STR(x(i,a))
         END DO
         WRITE(fid,'(A)')

         WRITE(fid,'(4X,A)',ADVANCE='NO') "R:"
         DO i=1, dof
            IF (i .EQ. 1) THEN
               WRITE(fid,'(1X,1pE25.18)') R(i,a)
            ELSE
               WRITE(fid,'(7X,1pE25.18)') R(i,a)
            END IF
         END DO

         DO i=rowPtr(a), rowPtr(a+1)-1
            b = colPtr(i)
            WRITE(fid,'(4X,A)') REPEAT('-',36)
            WRITE(fid,'(4X,A)') "Col: "//STR(b)//" gcol: "//STR(ltg(b))
            WRITE(fid,'(4X,A)',ADVANCE='NO') "R: "
            DO j=1, dof
               WRITE(fid,'(A)',ADVANCE='NO') " "//STR(R(j,b))
            END DO
            WRITE(fid,'(A)')

            DO k=1, dof
               WRITE(fid,'(4X,A)',ADVANCE='NO')
               DO j=1, dof
                  WRITE(fid,'(A)',ADVANCE='NO') " "//
     2               STR(Val((k-1)*dof+j,i))
               END DO
               WRITE(fid,'(A)')
            END DO
         END DO
      END DO
      CLOSE(fid)

      RETURN
      END SUBROUTINE PDEBUGVALR
!--------------------------------------------------------------------
      SUBROUTINE DEBUGVALR(K, U)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: K(dof*dof,lhs%nnz), U(dof,tnNo)

      INTEGER(KIND=IKIND) a, b, i, Ac, iM, nNo, gnNo, fid
      REAL(KIND=RKIND) :: KU(dof,tnNo)
      CHARACTER(LEN=stdL) fName

      REAL(KIND=RKIND), ALLOCATABLE :: lX(:,:), gX(:,:), lR(:,:),
     2   gR(:,:), lU(:,:), gU(:,:)

      KU(:,:) = 0._RKIND
      SELECT CASE (dof)
      CASE (1)
         DO a=1, tnNo
            DO i=rowPtr(a), rowPtr(a+1)-1
               b = colPtr(i)
               KU(1,a) = KU(1,a) + K(1 ,i)*U(1,b) + K(2 ,i)*U(2,b) +
     2                             K(3 ,i)*U(3,b) + K(4 ,i)*U(4,b)
            END DO
         END DO
      CASE (2)
         DO a=1, tnNo
            DO i=rowPtr(a), rowPtr(a+1)-1
               b = colPtr(i)
               KU(1,a) = KU(1,a) + K(1 ,i)*U(1,b) + K(2 ,i)*U(2,b) +
     2                             K(3 ,i)*U(3,b) + K(4 ,i)*U(4,b)
               KU(2,a) = KU(2,a) + K(5 ,i)*U(1,b) + K(6 ,i)*U(2,b) +
     2                             K(7 ,i)*U(3,b) + K(8 ,i)*U(4,b)
            END DO
         END DO
      CASE (3)
         DO a=1, tnNo
            DO i=rowPtr(a), rowPtr(a+1)-1
               b = colPtr(i)
               KU(1,a) = KU(1,a) + K(1 ,i)*U(1,b) + K(2 ,i)*U(2,b) +
     2                             K(3 ,i)*U(3,b) + K(4 ,i)*U(4,b)
               KU(2,a) = KU(2,a) + K(5 ,i)*U(1,b) + K(6 ,i)*U(2,b) +
     2                             K(7 ,i)*U(3,b) + K(8 ,i)*U(4,b)
               KU(3,a) = KU(3,a) + K(9 ,i)*U(1,b) + K(10,i)*U(2,b) +
     2                             K(11,i)*U(3,b) + K(12,i)*U(4,b)
            END DO
         END DO
      CASE (4)
         DO a=1, tnNo
            DO i=rowPtr(a), rowPtr(a+1)-1
               b = colPtr(i)
               KU(1,a) = KU(1,a) + K(1 ,i)*U(1,b) + K(2 ,i)*U(2,b) +
     2                             K(3 ,i)*U(3,b) + K(4 ,i)*U(4,b)
               KU(2,a) = KU(2,a) + K(5 ,i)*U(1,b) + K(6 ,i)*U(2,b) +
     2                             K(7 ,i)*U(3,b) + K(8 ,i)*U(4,b)
               KU(3,a) = KU(3,a) + K(9 ,i)*U(1,b) + K(10,i)*U(2,b) +
     2                             K(11,i)*U(3,b) + K(12,i)*U(4,b)
               KU(4,a) = KU(4,a) + K(13,i)*U(1,b) + K(14,i)*U(2,b) +
     2                             K(15,i)*U(3,b) + K(16,i)*U(4,b)
            END DO
         END DO
      END SELECT

      CALL COMMU(KU)

      fid = 1256
      i   = eq(cEq)%itr
      WRITE(fName,'(A)') TRIM(appPath)//"dbgValR_"//TRIM(eq(cEq)%sym)//
     2   "_"//STR(cTS)//"_"//STR(i)
      IF (cm%mas()) OPEN(fid,FILE=TRIM(fName))

      DO iM=1, nMsh
         nNo  = msh(iM)%nNo
         gnNo = msh(iM)%gnNo
         ALLOCATE(lX(nsd,nNo), lR(dof,nNo), lU(dof,nNo))
         IF (cm%mas()) THEN
            ALLOCATE(gX(nsd,gnNo), gR(dof,gnNo), gU(dof,gnNo))
         ELSE
            ALLOCATE(gX(0,0), gR(0,0), gU(0,0))
         END IF
         DO a=1, nNo
            Ac = msh(iM)%gN(a)
            lX(:,a) = x (:,Ac)
            lR(:,a) = R (:,Ac)
            lU(:,a) = KU(:,Ac)
         END DO

         gX = GLOBAL(msh(iM), lX)
         gR = GLOBAL(msh(iM), lR)
         gU = GLOBAL(msh(iM), lU)

         IF (cm%mas()) THEN
            WRITE(fid,'(A)') REPEAT('=',36)
            WRITE(fid,'(A)') "Mesh: <"//TRIM(msh(iM)%name)//">"
            DO a=1, gnNo
               WRITE(fid,'(A)',ADVANCE='NO') STR(a)
               DO i=1, nsd
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gX(i,a))
               END DO
               DO i=1, dof
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gR(i,a))
               END DO
               DO i=1, dof
                  WRITE(fid,'(A)',ADVANCE='NO') " "//STR(gU(i,a))
               END DO
               WRITE(fid,'(A)')
            END DO
         END IF
         DEALLOCATE(lX, gX, lR, lU, gR, gU)
      END DO

      IF (cm%mas()) CLOSE(fid)

      iM = cm%reduce(iM)

      RETURN
      END SUBROUTINE DEBUGVALR
!####################################################################

max_stars_repo_name: ktbolt/svFSI

max_issues_repo_name: ktbolt/svFSI

max_forks_repo_name: ktbolt/svFSI


text:
      function commit_particles_src()
      include 'common.h'

      call sort2(nt, r, iname)

      end function

      function recommit_particles_src()
      include 'common.h'

      call sort2(nt, r, iname)

      end function

      function get_number_of_particles_src(n_)
      include 'common.h'
      integer get_number_of_particles_src
      integer n_

      n_ = NT
      get_number_of_particles_src = 0
      end function

      function get_time_src(time_)
      include 'common.h'
      integer get_time_src
      double precision timet_
      time_ = time
      get_time_src = 0
      end function

      function get_timet_src(timet_)
      include 'common.h'
      integer get_timet_src
      double precision timet_
      timet_ = timet
      get_timet_src = 0
      end function

      function get_crossing_time_src(tcr_)
      include 'common.h'
      integer get_crossing_time_src
      double precision tcr_
      tcr_ = tcr
      get_timet_src = 0
      end function

      function get_state_src(id, mass_, r_, vr_, vt_)
      include 'common.h'
      integer get_state_src
      integer id
      double precision mass_, r_, vr_, vt_

      mass_ = BODY(id)

      r_ = R(id)
      vr_ = VR(id)
      vt_ = VT(id)
      get_state_src = 0
      end function get_state_src

      function set_state_src(id, mass_, r_, vr_, vt_)
      include 'common.h'
      integer set_state_src
      integer id
      double precision mass_, r_, vr_, vt_

      BODY(id) = mass_

      R(id) = r_
      VR(id) = vr_
      VT(id) = vt_
      set_state_src = 0
      end function set_state_src

      function get_total_kinetic_energy_src(T)
      include 'common.h'
      integer get_total_kinetic_energy
      double precision T
      T = ZKIN
      get_total_kinetic_energy = 0
      end function

      function get_total_potential_energy_src(V)
      include 'common.h'
      integer get_total_potential_energy_src
      double precision V
      V = -1
      get_total_potential_energy_src = 0
      end function

      function set_flagr(temp_, N)
      include 'common.h'
      integer, intent(in)::N
      double precision temp_(N)
      integer set_flagr, i
      do i = 1, N
          flagr(i) = temp_(i)
      end do
      set_flagr = 0
      end function

      function get_flagr(id, flagr_)
      include 'common.h'
      integer id
      double precision flagr_
      integer get_flagr
      flagr_ = flagr(id)
      get_flagr = 0
      end function

      function get_nlagra(temp_)
      include 'common.h'
      integer get_nlagra
      integer temp_
      temp_ = nlagra
      get_nlagra = 0
      end function

      function test_sort(n, aa, bb)
      integer n
      real*8 aa
      integer bb
      integer test_sort
      dimension aa(n), bb(n)

      call sort2(n, aa, bb)

      test_sort = 0
      end function

      function get_rc_src(rc_)
      include 'common.h'
      integer get_rc_src
      double precision rc_

      rc_ = rc

      get_rc_src = 0
      end function

      function get_nc_src(nc_)
      include 'common.h'
      integer get_rc_src
      integer nc_

      nc_ = nc

      get_nc_src = 0
      end function

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
*dk,recon
      subroutine recon(imsgin,xmsgin,cmsgin,msgtype,nwds,ierror)
C
C
C #####################################################################
C
C     PURPOSE -
C
C        This routine call the correct grid reconnection routine
C           depending on the type and dimensionality of the CMO.
C
C     INPUT ARGUMENTS -
C
C        imsgin()  - Integer array of command input tokens
C        xmsgin()  - Real array of command input tokens
C        cmsgin()  - Character array of command input tokens
C        msgtype() - Integer array of command input token types
C        nwds      - Number of command input tokens
C
C     OUTPUT ARGUMENTS -
C
C        ierror - Error Return Code (==0 ==> OK, <>0 ==> Error)
C
C     CHANGE HISTORY -
C
C        $Log: recon.f,v $
C        Revision 2.00  2007/11/09 20:04:00  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.9   24 Apr 2003 17:12:48   dcg
CPVCS    add check for ivornoi==5 which means recon has been disabled
CPVCS
CPVCS       Rev 1.8   15 May 2001 15:07:42   kuprat
CPVCS    We now compute TOLDAMAGE if user has not provided it, and
CPVCS    we pass it on to RECONLOOP2D and RECONLOOP3D.  We accept
CPVCS    the key word CHECKAXY and, if provided, pass a corresponding
CPVCS    flag on to RECONLOOP2D so that new triangles formed will
CPVCS    have an xy-projected area greater than EPSILONA.
CPVCS
CPVCS       Rev 1.7   12 Jun 2000 08:37:14   jan
CPVCS    changed line 230 from idum= ' ' to idum=0
CPVCS
CPVCS       Rev 1.6   Thu Apr 06 13:44:00 2000   dcg
CPVCS    replace get_info_i call
CPVCS
CPVCS       Rev 1.5   Wed Feb 03 15:24:02 1999   dcg
CPVCS    remove calls to fluxing routines and associated memory.
CPVCS
CPVCS       Rev 1.4   Fri Apr 17 15:35:28 1998   dcg
CPVCS    separate recon routines
CPVCS
CPVCS       Rev 1.61   Thu Apr 16 12:44:12 1998   dcg
CPVCS    add code to do 4to4 interface and volume flips if
CPVCS    ivoronoi = 2 or -2
CPVCS
CPVCS       Rev 1.59   Fri Feb 13 13:37:58 1998   kuprat
CPVCS    Deleted a "control C" that somehow got into the file.
CPVCS
CPVCS       Rev 1.58   Fri Feb 13 13:34:46 1998   kuprat
CPVCS    Commented out section writing to file called reconlist.
CPVCS
CPVCS       Rev 1.57   Mon Dec 01 16:41:22 1997   dcg
CPVCS    remove ipid
CPVCS    declare ipiholes as a pointer
CPVCS
CPVCS       Rev 1.56   Thu Oct 16 14:31:58 1997   kuprat
CPVCS    RECON format is now
CPVCS    recon/{0,1}/toldamage
CPVCS    2to2r flips causing boundary 'damage' exceeding TOLDAMAGE
CPVCS    are now prohibited.
CPVCS
CPVCS       Rev 1.55   Thu Oct 09 09:52:24 1997   dcg
CPVCS    fix memory management calls for merglst and merglst2
CPVCS
CPVCS       Rev 1.54   Wed Oct 08 16:55:58 1997   dcg
CPVCS    add error flag argument to one of the flip calls
CPVCS
CPVCS       Rev 1.53   Wed Oct 08 14:37:24 1997   dcg
CPVCS    fix typo on format statement
CPVCS
CPVCS       Rev 1.52   Wed Oct 08 13:55:20 1997   dcg
CPVCS    replace print statements
CPVCS
CPVCS       Rev 1.51   Wed Oct 01 14:40:14 1997   dcg
CPVCS    use mmfindbk with merge temp storage
CPVCS
CPVCS       Rev 1.50   Mon Sep 29 15:01:42 1997   dcg
CPVCS    call cmo_newlen at end of reconloop3d not at end of recon2
CPVCS
CPVCS       Rev 1.49   Tue Sep 16 16:04:58 1997   dcg
CPVCS    release merglst block
CPVCS
CPVCS       Rev 1.48   Wed Jul 09 09:12:34 1997   dcg
CPVCS    recheck insphere test at each iteration
CPVCS    flips from previous iterations may change results of insphere test
CPVCS    improve diagnostics if idebug is set
CPVCS
CPVCS       Rev 1.47   Fri Jun 06 11:27:00 1997   dcg
CPVCS    fix memory problems with kfix and xfix
CPVCS
CPVCS       Rev 1.46   Mon Jun 02 16:02:56 1997   dcg
CPVCS    change test to increment kfix allocation to account
CPVCS    for possible increase in ntets of 100 during next
CPVCS    iteration
CPVCS
CPVCS       Rev 1.45   Fri May 16 10:31:04 1997   kuprat
CPVCS    Added guard on 'N5' calculation.
CPVCS
CPVCS       Rev 1.44   Wed May 14 12:44:26 1997   dcg
CPVCS    fix memory problem with kfix overwritten if
CPVCS    3 to 2 flips increase ntets
CPVCS    fix output line to make more sense
CPVCS
CPVCS       Rev 1.43   Mon May 05 17:18:50 1997   dcg
CPVCS    resolve temporary storage conflicts between
CPVCS    recon and merge
CPVCS
CPVCS       Rev 1.42   Thu Apr 17 16:04:46 1997   dcg
CPVCS    set nelements by calling cmo_set_info before
CPVCS    call to tettestd
CPVCS
CPVCS       Rev 1.0   Thu Apr 17 16:03:32 1997   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.41   Mon Apr 14 16:58:22 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.40   Wed Apr 02 11:38:02 1997   dcg
CPVCS    change default 3d recon to recon/0
CPVCS
CPVCS       Rev 1.39   Tue Feb 04 03:30:30 1997   het
CPVCS    Add the neibor.h include file to reconloop3d.
CPVCS
CPVCS       Rev 1.38   Fri Jan 24 13:39:06 1997   het
CPVCS    Change the matrecon()=1 values and put unreconnected connections
CPVCS       into a file called reconlist.
CPVCS
CPVCS       Rev 1.37   Tue Dec 03 13:18:10 1996   het
CPVCS    Write out a list of unreconnect faces to a file.
CPVCS
CPVCS       Rev 1.36   Tue Nov 12 09:54:40 1996   dcg
CPVCS    test for ivoronoi =2 to avoid loop to fix
CPVCS    negative coupling coefficients
CPVCS
CPVCS       Rev 1.35   Thu Oct 17 16:12:14 1996   dcg
CPVCS
CPVCS
CPVCS
CPVCS
CPVCS    add ivoronoi=-2 option to set igeom_gwt to 1
CPVCS
CPVCS       Rev 1.34   Mon Jul 15 12:37:46 1996   dcg
CPVCS    fix call to mmfindbk
CPVCS
CPVCS       Rev 1.33   Thu Jun 27 15:26:56 1996   dcg
CPVCS    add constrainv command to create xcontab table as needed
CPVCS
CPVCS       Rev 1.32   Wed May 22 10:16:30 1996   dcg
CPVCS    get nconbnd from mesh object
CPVCS
CPVCS       Rev 1.31   Tue May 07 17:23:52 1996   het
CPVCS    Add the GWT reconnection algorithm.
CPVCS
CPVCS       Rev 1.30   Tue Apr 16 14:36:38 1996   dcg
CPVCS    try2to2b.f try2to4r.f try4to4x.f trymtonr.f
CPVCS    replace pointer ipitets with array itets
CPVCS
CPVCS       Rev 1.29   Tue Apr 02 02:29:28 1996   het
CPVCS    Fix an error the jtet1(jtet1)=something.
CPVCS
CPVCS       Rev 1.28   Sat Mar 09 10:29:00 1996   dcg
CPVCS    set t1,t2 to zero note timing not implemented
CPVCS
CPVCS       Rev 1.27   Tue Mar 05 12:37:46 1996   het
CPVCS    Move the (cmo)cmo_get_info calles from the driver into this routine.
CPVCS
CPVCS       Rev 1.26   Tue Feb 13 14:56:30 1996   het
CPVCS    Correct an error with setting interface mbndry to the external mbndry flag.
CPVCS
CPVCS       Rev 1.25   Fri Feb 02 14:23:20 1996   dcg
CPVCS    remove references to explicit vector attributes (u,w,v,e,r,pic)
CPVCS
CPVCS       Rev 1.24   Fri Dec 22 14:11:08 1995   het
CPVCS    Extract the iopt2to2 variable from the dictionary.
CPVCS
CPVCS       Rev 1.23   11/17/95 15:55:04   dcg
CPVCS    fix set_info_i calls
CPVCS
CPVCS       Rev 1.22   11/16/95 15:22:12   dcg
CPVCS    replace character literals in calls
CPVCS
CPVCS       Rev 1.21   11/07/95 17:23:36   dcg
CPVCS    change flag to 2 in mmgetblk calls
CPVCS
CPVCS       Rev 1.20   11/07/95 11:28:50   het
CPVCS    Modify the 2D triangle refinement algorithms.
CPVCS
CPVCS       Rev 1.19   10/22/95 13:18:24   het
CPVCS    Correct a memory management error
CPVCS
CPVCS       Rev 1.18   10/20/95 15:49:30   dcg
CPVCS    disable debug print
CPVCS
CPVCS       Rev 1.17   09/29/95 09:13:30   het
CPVCS    Put in added attributes inheritance
CPVCS
CPVCS       Rev 1.16   08/30/95 21:09:12   het
CPVCS    Put cmo table data into the cmoatt storage block
CPVCS
CPVCS       Rev 1.15   08/29/95 12:02:52   dcg
CPVCS    set length for names to 40 characters
CPVCS
CPVCS       Rev 1.14   08/16/95 00:23:26   het
CPVCS    Correct errors in recon and refine calles
CPVCS
CPVCS       Rev 1.13   08/15/95 18:23:50   het
CPVCS    Cleanup code and correct errors
C
c
c.... Format:
c....         RECON/[0,1]/[toldamage]/[CHECKAXY]
c....
c.... 1=Reconnection WITH edge refinement if
c....   necessary to eliminate negative coupling coefficients.
c.... 0=We only do reconnections, leaving any negative
c....   coupling coeffients that may still exist afterwards.
c.... TOLDAMAGE= maximum normal deformation allowed by reconnection.
c.... CHECKAXY:  If .true., then for the case of 2-D triangular
c....            meshes, we check xy- projected areas are positive
c....            and larger than EPSILONA.
C ######################################################################
C
      implicit none
C
      include "local_element.h"
      include "chydro.h"
      include "neibor.h"
      include "consts.h"

C arguments
      integer nwds, imsgin(nwds), msgtype(nwds)
      real*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)

C variables
      pointer (ipimt1, imt1)
      pointer (ipitp1, itp1)
      integer imt1(*), itp1(*)

      pointer (ipitet, itet )
      integer itet(*)

      pointer (ipitlist, itlist)
      integer itlist(*)

      pointer (iptable2, table2 )
      integer table2(*)

      pointer (ipxcontab, xcontab )
      real*8 xcontab(*)

      integer ierror,idum,ierrwrt,ijob,ilen,itype,icscode,nen,nef,nsd
     &   ,nelements,length,icmotype,ier,nelements_save,npoints,l,i
     &   ,nconbnd,ntets,it,lenfix,ntlist,i1,ics,lenout,icharlnf

      real*8 toldamage,xmin,xmax,ymin,ymax,zmin,zmax,epsilona,epsilonv

      logical lcheckaxy

      character*132 logmess
      character*32 isubname,contblnm
      character*32 cmo
      character*8 cdum
C
C#######################################################################
C BEGIN begin
C
      isubname="recon"
      contblnm='contable'
      cdum='-notset-'
      idum=0
C
      call cmo_get_name(cmo,ierror)
      if(ierror.ne.0) then
         write(logmess,9000)
 9000    format('No CMOs defined')
         call writloga('default',1,logmess,1,ierrwrt)
         goto 9999
      endif
      if(nwds.ge.2) then
         ijob=imsgin(2)
      else
         ijob=0
      endif

C     check for 0 elements
      call cmo_get_info('nelements',cmo,ntets,length,icmotype,ierror)
      if (ntets .le. 0) then
         write(logmess,'(a)')
     *   'WARNING Recon Early Exit: 0 elements'
         call writloga('default',1,logmess,1,ierror)
         ierror = -1
         goto 9999
      endif

C
c.... If provided, get damage tolerance from argument list.
c.... Otherwise, we get the dimensions of mesh object (using setsize and getsize)
c.... and we set TOLDAMAGE equal to .01 times the diameter of the
c.... mesh object.

      call setsize()
      call getsize(xmin,xmax,ymin,ymax,zmin,zmax,epsilona,epsilonv)
      if (msgtype(3).gt.0.and.cmsgin(3)(1:5).ne.'-def-') then
         call test_argument_type(1,2,3,imsgin,xmsgin,cmsgin,msgtype,
     *      nwds)
         toldamage=xmsgin(3)
      else
         toldamage=.01*sqrt((xmax-xmin)**2+(ymax-ymin)**2+
     &      (zmax-zmin)**2)
      endif

c.... If CHECKAXY is given, we make sure the projected areas of the
c.... triangles (in the case of 2-D flipping) are positive and
c.... greater than EPSILONA.

      lcheckaxy=.false.
      if (msgtype(4).gt.0.and.cmsgin(4)(1:5).ne.'-def-') then
         call test_argument_type(1,3,4,imsgin,xmsgin,cmsgin,msgtype,
     *      nwds)
         if (cmsgin(4)(1:8).eq.'checkaxy') then
            lcheckaxy=.true.
         endif
      endif

      call cmo_get_info('ivoronoi',cmo,
     *                ivoronoi,ilen,itype,icscode)
         if (icscode .ne. 0) call x3d_error(isubname,'cmo_get_info')
      if(ivoronoi.eq.5) go to 9999
      if(abs(ivoronoi).eq.2) then
         call mega_hessian()
         call mega_error()
      endif
      if(ivoronoi.eq.2.and.ijob.eq.1) ijob=0
C
      call cmo_get_info('nodes_per_element',cmo,nen,ilen,itype,icscode)
      call cmo_get_info('faces_per_element',cmo,nef,ilen,itype,icscode)
      call cmo_get_info('ndimensions_topo',cmo,nsd,ilen,itype,icscode)
C
      call cmo_get_info('nelements',cmo,nelements,length,icmotype,ier)
      nelements_save=nelements
C
      if(nen.eq.nelmnen(ifelmtri).and.nef.eq.nelmnef(ifelmtri)) then
         call reconloop2d(ijob,toldamage,lcheckaxy,epsilona)

      elseif(nen.eq.nelmnen(ifelmtet).and.nef.eq.nelmnef(ifelmtet).and.
     *       nsd.eq.3) then
C
         call cmo_get_name(cmo,ierror)
         call cmo_get_info('nnodes',cmo,npoints,l,i,ierror)
         call cmo_get_info('nconbnd',cmo,nconbnd,l,i,ierror)
         call cmo_get_info('nelements',cmo,ntets,l,i,ierror)
         call cmo_get_info('imt1',cmo,ipimt1,l,i,ierror)
         call cmo_get_info('itp1',cmo,ipitp1,l,i,ierror)
         call cmo_get_info('itet',cmo,ipitet,l,i,ierror)
C  set up temporary storage
C  set 1000 as min length to avoid single elements lacking memory
         l=max(l,4*ntets)
         length=max(1000,(l-1)/4+1)
         call mmgetblk("ifittet",isubname,ipifittet,length,2,icscode)
         do it=1,ntets
            ifittet(it)=0
         enddo
         length=max(1000,nint(1.7*4*length))
         call mmfindbk('kfix',isubname,ipkfix,lenfix,icscode)
         if(icscode.ne.0) then
            call mmgetblk('kfix',isubname,ipkfix,length,1,icscode)
         else
            if(lenfix.lt.length)
     *      call mmnewlen('kfix',isubname,ipkfix,length,icscode)
         endif
         call mmfindbk('xfix',isubname,ipxfix,lenfix,icscode)
         if(icscode.ne.0) then
            call mmgetblk('xfix',isubname,ipxfix,length,2,icscode)
         else
            if(lenfix.lt.length)
     *      call mmnewlen('xfix',isubname,ipxfix,length,icscode)
         endif
         do i=1,4*ntets
            kfix1(i)=0
         enddo
         ntlist=max(ntets,(l-1)/4+1)
         length=ntlist
         call mmgetblk("itlist",isubname,ipitlist,length,2,icscode)
         do it=1,ntets
            itlist(it)=it
         enddo
         do i=1,ntlist
            kfix1(i)=itlist(i)
         enddo
C
         length=6*npoints
         call mmgetblk('iinter',isubname,ipinter,length,2,icscode)
         ninter=0
         do i1=1,npoints
            if(itp1(i1).eq.ifitpcup) then
               ninter=ninter+1
               iinter(1,ninter)=i1
            endif
         enddo
C
         nmats=0
         do i1=1,npoints
            nmats=max(nmats,imt1(i1))
         enddo
         length=nmats
         call mmgetblk('matrecon',isubname,ipmatrecon,length,2,ics)
         do i=1,nmats
            matrecon(i)=1
         enddo
         call initnbr()
         if (nconbnd .gt. 0) then
            call mmfindbk('mtable2',contblnm,iptable2,lenout,ics)
            if(icscode.ne.0) call contable()
         endif
C
         call cmo_get_name(cmo,ierror)
C
C
         call reconloop3d(ijob,toldamage)
C
         call mmrelblk("itlist",isubname,ipitlist,icscode)
         call mmrelblk("ifittet",isubname,ipifittet,icscode)
         call mmrelblk("kfix",isubname,ipkfix,icscode)
         call mmrelblk("xfix",isubname,ipxfix,icscode)
         call mmrelblk('iinter',isubname,ipinter,icscode)
         call mmrelblk('matrecon',isubname,ipmatrecon,icscode)
         if (nconbnd .gt. 0) then
            call mmrelblk('mtable2',contblnm,iptable2,icscode)
         endif
      else
         write(logmess,9010) cmo(1:icharlnf(cmo))
 9010    format('Refine on this CMO type is not implemented: ',a)
         call writloga('default',1,logmess,1,ierrwrt)
         goto 9999
      endif
C
      goto 9999
 9999 continue
C
C      CLEAR XCONTAB IF IT EXISTS
      call mmfindbk ('xcontab',cmo,ipxcontab,lenout,ier)
      if(ier.eq.0) call dotaskx3d('constrainv/clear ; finish',ier)
c
c  clear flip memory
c
      call mflip(99,idum,cdum)
C
      write(logmess,'(a)') 'RECON: done.'
      call writloga('default',0,logmess,0,ier)

      return
      end

max_stars_repo_name: daniellivingston/LaGriT-CI-Test

max_issues_repo_name: daniellivingston/LaGriT-CI-Test

max_forks_repo_name: daniellivingston/LaGriT-CI-Test


text:
bitwidth/stimulus.cpp
bitwidth/display.cpp
bitwidth/bitwidth.cpp
bitwidth/main.cpp


text:
************************************************************************
*
*     Subroutine AMLSTAT                         Called by: CALIBR
*
*     compute AMLE statistics and residuals for the best model
*
*     local vars
*     ----------
*     RESSD    residual standard deviation
*     RESPPCC  residuals for PPCC calculation
*     RESSORT  copy of residuals for use in SORT
*     SUM      sum of YHAT
*     SUMSQ    sum of squares, YHAT
*     SUMSR    sum of squared residuals
*     TAKEN    variable used to assign Z Scores
*
************************************************************************
      SUBROUTINE AMLSTAT(CENSFLAG,IPPCCERR,LDOUT,NOBSC,NOBSCI,NPAR,
     &                   PARAML,PPCC,PLEV,PVAL,RESID,RSQ,SCORR,XLCAL,
     &                   XLIKE,YD,YHAT,YLCAL,YLCAL2,Z)
*
*     subroutine arguments
*
      LOGICAL CENSFLAG(*)
      INTEGER*4 IPPCCERR,LDOUT,NOBSC,NOBSCI,NPAR
      DOUBLE PRECISION PPCC,PLEV,RSQ,SCORR,XLIKE
      DOUBLE PRECISION PARAML(*),PVAL(*),RESID(*),YD(*),YHAT(*),
     &                 YLCAL(*),YLCAL2(*),Z(*)
      DOUBLE PRECISION XLCAL(NOBSCI,*)
*
*     local vars
*
      INTEGER*4 I,K
      DOUBLE PRECISION A,RESSD,SUM,SUMSR,SUMSQ
      DOUBLE PRECISION RESPPCC(NOBSCI),RESSORT(NOBSCI),ZWORK(NOBSCI)
      LOGICAL TAKEN(NOBSCI)
*
*     function declarations
*
      DOUBLE PRECISION CORR,PRED,ZSCORE
*
*     calculate P-values (PVAL) for the AMLE regression coefficients
*
      CALL AMLPVAL(CENSFLAG,LDOUT,NOBSC,NOBSCI,NPAR,PVAL,XLCAL,XLIKE,
     &             YD,YLCAL)
*
*     calculate residuals (observed load - predicted), SUM, and
*     SUMSQ.  For censored observations, use the expected value of
*     observations rather than the observations (5/7/92 Tim Cohn).
*
*     determine RESPPCC, the residual for the PPCC calculation; the
*     detection limit, rather than the expected value, is used to
*     calculate the residual when an observation is censored.
*
*     (Preserve the original YLCAL for use by LADREG and JKNIFE by
*     introducing YLCAL2).
*
      RESSD = SQRT(PARAML(NPAR+1))
      SUM = 0.D0
      SUMSQ = 0.D0
      DO 10 I=1,NOBSC
         YLCAL2(I) = YLCAL(I)
         YHAT(I) = PRED(NPAR,NOBSCI,I,XLCAL,PARAML)
         IF (CENSFLAG(I)) THEN
            CALL TACIT_CALC((YD(I)-YHAT(I))/RESSD,A)
            YLCAL2(I) = YHAT(I) - RESSD*A
         ENDIF
         RESID(I) = YLCAL2(I) - YHAT(I)
         RESPPCC(I) = YLCAL(I) - YHAT(I)
         SUM = SUM + YHAT(I)
         SUMSQ = SUMSQ + YHAT(I)**2
 10   CONTINUE
*
*     compute R^2 of the load regression and the serial correlation of
*     the residuals
*
      SUMSR = SUMSQ-SUM**2/NOBSC
      RSQ = 100.D0 * SUMSR/(((NOBSC-NPAR)*RESSD**2+SUMSR))
      SCORR = CORR(NOBSC-1,RESID(1),RESID(2))
*
*     make a copy of the residuals vector so that it may be sorted for
*     use in assigning Z-scores, determine Z-scores, and initialize
*     TAKEN
*
      DO 20 I=1,NOBSC
         RESSORT(I) = RESID(I)
         ZWORK(I) = ZSCORE(DBLE(I)/DBLE(NOBSC+1))
         TAKEN(I) = .FALSE.
 20   CONTINUE
*
*     sort residuals
*
      CALL SORT(RESSORT,NOBSC)
*
*     output residuals, Z-scores, predicted values, flow, and decimal
*     time.  Prior to output, assign Z-scores.
*
*     (Assignment of Z-scores is necessary as Z corresponds to RESSORT,
*     rather than the original vector RESID.  The variable TAKEN is
*     included to make sure that residuals with identical magnitudes
*     are given different Z-scores).
*
      DO 40 I=1,NOBSC
         DO 30 K=1,NOBSC
            IF (RESID(I) .EQ. RESSORT(K) .AND. (.NOT. TAKEN(K))) THEN
               TAKEN(K) = .TRUE.
               Z(I) = ZWORK(K)
               GOTO 40
            ENDIF
 30      CONTINUE
 40   CONTINUE
*
*     compute the probability plot correlation coefficient
*
      CALL PPCCTEST(CENSFLAG,NOBSC,RESPPCC,PPCC,PLEV,IPPCCERR)

      RETURN
      END

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


text:
! ********** HILFSROUTINEN **********

module fdschemes

  implicit none

contains

  ! ***** Forward Algorithmus *****
  subroutine forward(u0,u,a,dx,dt)

    ! Eingabe:
    ! u0 = Ausgangsvektor
    ! u = Ausgabevektor
    ! a = Geschwindigkeit (>0)
    ! dx = Ortsschrittweite
    ! dt = Zeitschrittweite

    double precision, intent(in)			:: a, dx, dt
    double precision, dimension(:), intent(out)		:: u
    double precision, dimension(:), intent(in)		:: u0
    integer						:: i,j
    double precision					:: k1

    ! Hilfsvariable
    k1 = a*dt/dx

    ! Berechne den Loesungsvektor fuer alle ausser dem linken x (parallel)
    forall(i=lbound(u0,1)+1:ubound(u0,1))
       u(i) = u0(i)-k1*(u0(i)-u0(i-1))
    end forall

    ! Berechne den Loesungswert fuer das linke x
    u(lbound(u0,1)) = u0(lbound(u0,1))+k1*u0(lbound(u0,1))

  end subroutine forward

  ! ***** generalisierter Upwind und Lax-Wendroff *****
  subroutine gul(u0, u, x, dt, dx, funcf, funcfstrich, alg)
	! Loest:
	! d/dt u + d/dx f(u) = 0
	! u(x,0) = u0(x)

	! Eingabeparameter:
	! u0            Startwert-Vektor (Zeilenvektor)
	! x		Ortsdiskretisierung (Zeilenvektor)
	! dt            Laenge des Zeitschritts
	! dx            Schrittweite des Ortsgitters
	! funcf         Funktionhandle auf die Funktion f(u)
	! funcfstrich   Funktionhandle auf die Funktion f'(u)
	! alg           waehlt den zu verwendenden Algorithmus
	!               0: generalisierter Upwind
	!               1: generalisierter Lax Wendroff

	double precision, intent(in)				:: dx, dt, alg
	double precision, dimension(:), intent(out)		:: u
	double precision, dimension(:), intent(in)		:: u0, x
	integer							:: i,j
	double precision					:: k1

	interface
	double precision function funcf(u,i,x)
		double precision, intent(in)	:: u, x
		integer, intent(in)		:: i
    	end function funcf
    	double precision function funcfstrich(u,i,x)
		double precision, intent(in)	:: u, x
		integer, intent(in)		:: i
    	end function funcfstrich
	end interface

	! Loesungsvektor: Randbedingung
	u(lbound(u0,1)) = u0(lbound(u0,1))
	u(ubound(u0,1)) = u0(ubound(u0,1))

	! Hilfsvariablen
	k1 = dt/dx

	! Berechne alle inneren Werte

	gulloop: do i = lbound(u0,1)+1,ubound(u0,1)-1
		u(i) = u0(i) - k1*(fl(i)-fl(i-1))
	end do gulloop

  contains

  ! Hilfsfunktionen fuer gul
  ! flux-function f_{i+1/2}  =f^L fuer alpha = 0, f^H fuer alpha = 1
    double precision function fl(i)
	double precision	:: k2, tempa
	integer, intent(in) 	:: i
        tempa = approxa(i);
        k2 = tempa*k1;    ! = nu
        fl = (funcf(u0(i+1),i+1,x(i+1))+funcf(u0(i),i,x(i)))/2.0d0-abs(tempa)/2.0d0*(u0(i+1)-u0(i))&
		+alg*abs(tempa)/2.0d0*(1-abs(k2))*(u0(i+1)-u0(i))
    end function fl
  ! berechnet a_{i+1/2}
    double precision function approxa(i)
	integer, intent(in)	:: i
        if (abs(u0(i+1)-u0(i))< 0.001d0) then
            approxa = funcfstrich(u0(i),i,x(i))
        else
            approxa = (funcf(u0(i+1),i+1,x(i+1))-funcf(u0(i),i,x(i)))/(u0(i+1)-u0(i))
        end if
    end function approxa

 end subroutine gul

  ! ***** generalisierter Upwind und Lax-Wendroff *****
  subroutine gtvd(u0, u, x, dt, dx, funcf, funcfstrich, limiter)
	! Loest:
	! d/dt u + d/dx f(u) = 0
	! u(x,0) = u0(x)

	! Eingabeparameter:
	! u0            Startwert-Vektor (Zeilenvektor)
	! x		Ortsdiskretisierung (Zeilenvektor)
	! dt            Laenge des Zeitschritts
	! dx            Schrittweite des Ortsgitters
	! funcf         Funktionhandle auf die Funktion f(u)
	! funcfstrich   Funktionhandle auf die Funktion f'(u)
	! limiter	waehlt den zu verwendenden Limiter (1=minmod, 2=mc, 3=vl, 4=superbee)

	double precision, intent(in)				:: dx, dt
	double precision, dimension(:), intent(out)		:: u
	double precision, dimension(:), intent(in)		:: u0, x
	integer, intent(in)					:: limiter
	integer							:: i,j
	double precision					:: k1

	interface
	double precision function funcf(u,i,x)
		double precision, intent(in)	:: u, x
		integer, intent(in)		:: i
    	end function funcf
    	double precision function funcfstrich(u,i,x)
		double precision, intent(in)	:: u, x
		integer, intent(in)		:: i
    	end function funcfstrich
	end interface

	! Loesungsvektor: Randbedingungen
	u(lbound(u0,1)) = u0(lbound(u0,1))
	u(lbound(u0,1)+1) = u0(lbound(u0,1)+1)
	u(ubound(u0,1)) = u0(ubound(u0,1))
	u(ubound(u0,1)-1) = u0(ubound(u0,1)-1)

	! Hilfsvariablen
	k1 = dt/dx

	! Berechne alle inneren Werte
	tvdloop: do i = lbound(u0,1)+2,ubound(u0,1)-2
		u(i) = u0(i) - k1*(fl(i)-fl(i-1))
	end do tvdloop

  contains

  ! Hilfsfunktionen fuer gtvd
  ! flux-function f_{i+1/2}  =f^L fuer alpha = 0, f^H fuer alpha = 1 (alpha = erg. des limiters)
    double precision function fl(i)
	double precision	:: k2, tempa
	integer, intent(in) 	:: i
	integer			:: j

        tempa = approxa(i)
        k2 = tempa*k1    ! = nu

	j=i-int(sign(1d0,tempa))          !???????????????????????????????????????????????????????????????
	!if (sign(tempa,1.0d0)>0d0) then
	!j=i-1
	!else
	!j=1+1
	!end if

        fl = (funcf(u0(i+1),i+1,x(i+1))+funcf(u0(i),i,x(i)))/2.0&
		-abs(tempa)/2.0*(u0(i+1)-u0(i))&
		 + limiterfunc(u0(i+1)-u0(i),u0(j+1)-u0(j))*abs(tempa)/2.0*(1-abs(k2))
    end function fl
  ! berechnet a_{i+1/2}
    double precision function approxa(i)
	integer, intent(in)	:: i
        if (abs(u0(i+1)-u0(i))< 0.001d0) then
            approxa = funcfstrich(u0(i),i,x(i))
        else
            approxa = (funcf(u0(i+1),i+1,x(i+1))-funcf(u0(i),i,x(i)))/(u0(i+1)-u0(i))
        end if
    end function approxa
  ! limiter function
    double precision function limiterfunc(a,b)
	double precision, intent(in)		:: a,b
	double precision			:: temp5,na,nb
        temp5=(sign(1.0d0,a)+sign(1.0d0,b))/2.0d0   ! ???????????????????
        na=abs(a)
        nb=abs(b)
        if (temp5.ne.0d0) then
        select case(limiter)
            case (1)
                    limiterfunc = temp5*min(na,nb)
            case (2)
                    if (na+nb < 0.00001d0) then
                        limiterfunc=1
                    else
                        limiterfunc = temp5*2.0d0*na*nb/(na+nb)
                    end if
            case (3)
                    limiterfunc = temp5*min(2.0d0*na,(na+nb)/2.0d0,2.0d0*nb)
            case (4)
                    limiterfunc = temp5*max(min(2.0d0*na,nb),min(na,2.0d0*nb))
        end select
        else
            limiterfunc=0d0
        end if
    end function limiterfunc

 end subroutine gtvd

  ! ***** Benchmarkfunctions *****
    double precision function bench1funcf(u,i,x)
	double precision, intent(in)	:: u,x
	integer, intent(in)		:: i
        bench1funcf = x*u
    end function bench1funcf
    double precision function bench1funcfstrich(u,i,x)
	double precision, intent(in)	:: u,x
	integer, intent(in)		:: i
        bench1funcfstrich = x
    end function bench1funcfstrich
    double precision function bench2funcf(u,i,x)
	double precision, intent(in)	:: u,x
	integer, intent(in)		:: i
        bench2funcf = (0.5d0-x)*u
    end function bench2funcf
    double precision function bench2funcfstrich(u,i,x)
	double precision, intent(in)	:: u,x
	integer, intent(in)		:: i
        bench2funcfstrich = 0.5d0-x
    end function bench2funcfstrich
    double precision function bench3funcf(u,i,x)
	double precision, intent(in)	:: u,x
	integer, intent(in)		:: i
        bench3funcf = u*u/2d0
    end function bench3funcf
    double precision function bench3funcfstrich(u,i,x)
	double precision, intent(in)	:: u,x
	integer, intent(in)		:: i
        bench3funcfstrich = u
    end function bench3funcfstrich

end module fdschemes

! ********** HAUPTPROGRAMM **********

program finite_differences
  use fdschemes

  implicit none

  integer, parameter				:: n = 101		! Anzahl der x-Gitterpunkte
  double precision, parameter			:: tfinal = 0.4		! Endzeit
  double precision				:: a, dx, dt, time
  double precision				:: time0, time1		! Zur Zeitmessung
  double precision, dimension(1:n)		:: u0 ,u ,x		! Startvektor, Loesungsvektor, Gittervektor
  integer					:: i			! Laufvariable

  ! Parameter initialisieren
  a = 1.0				! Welleneschwindigkeit (fuer die linear advection equation)
  dx = 1.0d0/(n-1)			! dx
  dt = 0.001d0				! dt

  ! Gittervektor initialisieren
  forall(i=1:n)
     x(i)=(i-1)*dx
  end forall

  ! Startwerte festlegen
  u0 = 0
  where (x>=0.1d0 .AND. x<=0.4d0)
     u0 = 1
  end where

  ! Datei oeffnen
  open(10,file='ausgabe.dat')

  ! Startvektor in Datei schreiben
  write(10,*) u0

  ! Zeitschritte ausfuehren und Ergebnis in Datei schreiben
  ! mit Zeitmessung
  call cpu_time(time0)
  time = 0.0
  timestepping: do
     time = time + dt
     if (time>tfinal) exit timestepping
     call gtvd(u0, u, x, dt, dx, bench3funcf, bench3funcfstrich, 4)
     u0 = u
     write(10,*) u
  end do timestepping
  call cpu_time(time1)

  ! Datei schliessen
  close(10)

  ! Zeitaufwand ausgeben
  print*,'Zeitaufwand: ',time1-time0,' Sekunden.'
end program finite_differences

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
Some things Im interested in:

University Airport
Cal Aggie Flying Farmers
Friends of the University Airport
VCOA Volvo Show and Swap Meet VCOAs Annual Davis Volvo Show and Swap Meet
Yolo Masonic Lodge

when is the next Yolo swap meet? Daubert

20130926 18:13:16 nbsp Hey, quick question  is the name of the group Friends of University Airport or Friends of the University Airport? Users/JabberWokky Evan JabberWokky Edwards

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module subspace_solv_sl_module

  use wf_module
  use scalapack_module, only: allocated_workarray_scalapack, UPLO, DESCA, DESCZ, &
  NP0,NQ0,NPX,NQX,MBSIZE,NPCOL
  use subspace_diag_variables, only: MB_diag,Hsub,Vsub
  use watch_module

  implicit none

  private
  public :: subspace_solv_sl

contains

  subroutine subspace_solv_sl(k,s)
    implicit none
    integer,intent(IN) :: k,s
    integer :: itmp(1),LWORK0,LRWORK0,LIWORK0,TRILWMIN,ierr,MB
    integer,save :: LWORK=0,LRWORK=0,LIWORK=0
    integer,allocatable :: iwork(:)
    real(8) :: rtmp(1)
    real(8),allocatable :: rwork(:)
    complex(8) :: ctmp(1)
    complex(8),allocatable :: zwork(:)
    character(8) :: idiag
    type(time) :: t
    logical :: disp_sw, ltmp

#ifndef _LAPACK_

    call write_border( 1, " subspace_solv_sl(start)" )
    call start_timer( t_out=t )

    MB = MB_diag

    ierr = 0

    call check_disp_switch( disp_sw, 0 )

#ifdef _DRSDFT_
    idiag = 'PDSYEVD'
#else
    idiag = 'PZHEEVD'
#endif

    select case( idiag )
    case( 'PZHEEVD' )

      if ( .not.allocated_workarray_scalapack() ) then
        call PZHEEVD('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1 &
                    ,DESCZ,ctmp,-1,rtmp,-1,itmp,-1,ierr)
        LWORK =nint(real(ctmp(1)))
        LRWORK=nint(rtmp(1))
        LIWORK=itmp(1)
        ! if ( disp_sw ) then
        !   write(*,'("(PZHEEVD) Work-array sizes (by query)   : LWORK,LRWORK,LIWORK=",3i8)') LWORK,LRWORK,LIWORK
        !   write(*,'("(PZHEEVD) Work-array sizes (in document): LWORK,LRWORK,LIWORK=",3i8)') &
        !     MB+(NP0+NQ0+MBSIZE)*MBSIZE, 1+9*MB+3*NPX*NQX, 7*MB+8*NPCOL+2
        ! end if
        LWORK =max( LWORK , MB+(NP0+NQ0+MBSIZE)*MBSIZE )
        LRWORK=max( LRWORK, (1+9*MB+3*NPX*NQX) )
        LIWORK=max( LIWORK, 7*MB+8*NPCOL+2 )
        ! if ( disp_sw ) then
        !   write(*,'("(PZHEEVD) Work-array sizes (allocate)   : LWORK,LRWORK,LIWORK=",3i8)') LWORK,LRWORK,LIWORK
        ! end if
        ltmp = allocated_workarray_scalapack( .true. )
      end if

      allocate( zwork(LWORK),rwork(LRWORK),iwork(LIWORK) )

      call PZHEEVD('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1 &
                  ,DESCZ,zwork,LWORK,rwork,LRWORK,iwork,LIWORK,ierr)

      deallocate( iwork,rwork,zwork )

    case('pzheev')

      if ( LWORK==0 ) then
        call pzheev('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1 &
                    ,DESCZ,ctmp,-1,rtmp,-1,ierr)
          LWORK =nint(real(ctmp(1)))
          LRWORK=nint(rtmp(1))
      end if

      LWORK =max(LWORK,(NP0+NQ0+MBSIZE)*MBSIZE+3*MB+MB**2)
      LRWORK=max(LRWORK,4*MB-2)

      allocate( zwork(LWORK),rwork(LRWORK) )

      call pzheev('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1 &
            ,DESCZ,zwork,LWORK,rwork,LRWORK,ierr)

      deallocate( rwork,zwork )

    case( 'PDSYEVD' )

      TRILWMIN = 3*MB + max( MBSIZE*(NPX+1),3*MBSIZE )
      LRWORK0  = max( 1+6*MB+2*NPX*NQX, TRILWMIN )
      LIWORK0  = 7*MB+8*NPCOL+2

      if ( LRWORK==0 ) then
        call PDSYEVD('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1 &
                     ,DESCZ,rtmp,-1,itmp,LIWORK,ierr)
        LRWORK=nint(rtmp(1))
        LRWORK=LRWORK*10
      end if
      LRWORK=max(LRWORK,LRWORK0*10)
      LIWORK=max(LIWORK,LIWORK0)

      allocate( rwork(LRWORK),iwork(LIWORK) )

      call PDSYEVD('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1 &
                  ,DESCZ,rwork,LRWORK,iwork,LIWORK,ierr)

      deallocate( iwork,rwork )

    case('pdsyev')

       LRWORK0=0

       if ( LRWORK==0 ) then
          call pdsyev('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1,DESCZ,rtmp,-1,ierr)
          LRWORK=nint(rtmp(1))
       end if
       LRWORK=max(LRWORK,LRWORK0)

       allocate( rwork(LRWORK) )

       call pdsyev('V',UPLO,MB,Hsub,1,1,DESCA,esp(1,k,s),Vsub,1,1,DESCZ,rwork,LRWORK,ierr)

       deallocate( rwork )

    case default

       write(*,*) "idiag=",idiag
       write(*,*) "This routine is not available!"
       stop

    end select

    if ( ierr /= 0 ) then
       write(*,*) "ierr,idiag=",ierr,idiag
       stop
    end if

    call result_timer( "solv_sl", t )
    call write_border( 1, " subspace_solv_sl(end)" )

    return

#endif

  end subroutine subspace_solv_sl

end module subspace_solv_sl_module


!===-- module/iso_fortran_env.f90 ------------------------------------------===!
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
!===------------------------------------------------------------------------===!

! See Fortran 2018, clause 16.10.2
! TODO: These are placeholder values so that some tests can be run.

include '../include/flang/Runtime/magic-numbers.h' ! for IOSTAT= error/end code values

module iso_fortran_env

  use __Fortran_builtins, only: &
    event_type => __builtin_event_type, &
    lock_type => __builtin_lock_type, &
    team_type => __builtin_team_type

  implicit none

  integer, parameter :: atomic_int_kind = selected_int_kind(18)
  integer, parameter :: atomic_logical_kind = atomic_int_kind

  ! TODO: Use PACK([x],test) in place of the array constructor idiom
  ! [(x, integer::j=1,COUNT([test]))] below once PACK() can be folded.

  integer, parameter, private :: &
    selectedASCII = selected_char_kind('ASCII'), &
    selectedUCS_2 = selected_char_kind('UCS-2'), &
    selectedUnicode = selected_char_kind('ISO_10646')
  integer, parameter :: character_kinds(*) = [ &
    [(selectedASCII, integer :: j=1, count([selectedASCII >= 0]))], &
    [(selectedUCS_2, integer :: j=1, count([selectedUCS_2 >= 0]))], &
    [(selectedUnicode, integer :: j=1, count([selectedUnicode >= 0]))]]

  integer, parameter, private :: &
    selectedInt8 = selected_int_kind(2), &
    selectedInt16 = selected_int_kind(4), &
    selectedInt32 = selected_int_kind(9), &
    selectedInt64 = selected_int_kind(18),&
    selectedInt128 = selected_int_kind(38), &
    safeInt8 = merge(selectedInt8, selected_int_kind(0), &
                     selectedInt8 >= 0), &
    safeInt16 = merge(selectedInt16, selected_int_kind(0), &
                      selectedInt16 >= 0), &
    safeInt32 = merge(selectedInt32, selected_int_kind(0), &
                      selectedInt32 >= 0), &
    safeInt64 = merge(selectedInt64, selected_int_kind(0), &
                      selectedInt64 >= 0), &
    safeInt128 = merge(selectedInt128, selected_int_kind(0), &
                       selectedInt128 >= 0)
  integer, parameter :: &
    int8 = merge(selectedInt8, merge(-2, -1, selectedInt8 >= 0), &
                 digits(int(0,kind=safeInt8)) == 7), &
    int16 = merge(selectedInt16, merge(-2, -1, selectedInt16 >= 0), &
                  digits(int(0,kind=safeInt16)) == 15), &
    int32 = merge(selectedInt32, merge(-2, -1, selectedInt32 >= 0), &
                  digits(int(0,kind=safeInt32)) == 31), &
    int64 = merge(selectedInt64, merge(-2, -1, selectedInt64 >= 0), &
                  digits(int(0,kind=safeInt64)) == 63), &
    int128 = merge(selectedInt128, merge(-2, -1, selectedInt128 >= 0), &
                   digits(int(0,kind=safeInt128)) == 127)

  integer, parameter :: integer_kinds(*) = [ &
    selected_int_kind(0), &
    ((selected_int_kind(k), &
      integer :: j=1, count([selected_int_kind(k) >= 0 .and. &
                             selected_int_kind(k) /= &
                               selected_int_kind(k-1)])), &
     integer :: k=1, 39)]

  integer, parameter :: &
    logical8 = int8, logical16 = int16, logical32 = int32, logical64 = int64
  integer, parameter :: logical_kinds(*) = [ &
    [(logical8, integer :: j=1, count([logical8 >= 0]))], &
    [(logical16, integer :: j=1, count([logical16 >= 0]))], &
    [(logical32, integer :: j=1, count([logical32 >= 0]))], &
    [(logical64, integer :: j=1, count([logical64 >= 0]))]]

  integer, parameter, private :: &
    selectedReal16 = selected_real_kind(3, 4), &      ! IEEE half
    selectedBfloat16 = selected_real_kind(2, 37), &   ! truncated IEEE single
    selectedReal32 = selected_real_kind(6, 37), &     ! IEEE single
    selectedReal64 = selected_real_kind(15, 307), &   ! IEEE double
    selectedReal80 = selected_real_kind(18, 4931), &  ! 80x87 extended
    selectedReal64x2 = selected_real_kind(31, 307), & ! "double-double"
    selectedReal128 = selected_real_kind(33, 4931), & ! IEEE quad
    safeReal16 = merge(selectedReal16, selected_real_kind(0,0), &
                       selectedReal16 >= 0), &
    safeBfloat16 = merge(selectedBfloat16, selected_real_kind(0,0), &
                         selectedBfloat16 >= 0), &
    safeReal32 = merge(selectedReal32, selected_real_kind(0,0), &
                       selectedReal32 >= 0), &
    safeReal64 = merge(selectedReal64, selected_real_kind(0,0), &
                       selectedReal64 >= 0), &
    safeReal80 = merge(selectedReal80, selected_real_kind(0,0), &
                       selectedReal80 >= 0), &
    safeReal64x2 = merge(selectedReal64x2, selected_real_kind(0,0), &
                         selectedReal64x2 >= 0), &
    safeReal128 = merge(selectedReal128, selected_real_kind(0,0), &
                        selectedReal128 >= 0)
  integer, parameter :: &
    real16 = merge(selectedReal16, merge(-2, -1, selectedReal16 >= 0), &
                   digits(real(0,kind=safeReal16)) == 11), &
    bfloat16 = merge(selectedBfloat16, merge(-2, -1, selectedBfloat16 >= 0), &
                     digits(real(0,kind=safeBfloat16)) == 8), &
    real32 = merge(selectedReal32, merge(-2, -1, selectedReal32 >= 0), &
                   digits(real(0,kind=safeReal32)) == 24), &
    real64 = merge(selectedReal64, merge(-2, -1, selectedReal64 >= 0), &
                   digits(real(0,kind=safeReal64)) == 53), &
    real80 = merge(selectedReal80, merge(-2, -1, selectedReal80 >= 0), &
                   digits(real(0,kind=safeReal80)) == 64), &
    real64x2 = merge(selectedReal64x2, merge(-2, -1, selectedReal64x2 >= 0), &
                     digits(real(0,kind=safeReal64x2)) == 106), &
    real128 = merge(selectedReal128, merge(-2, -1, selectedReal128 >= 0), &
                    digits(real(0,kind=safeReal128)) == 113)

  integer, parameter :: real_kinds(*) = [ &
    [(real16, integer :: j=1, count([real16 >= 0]))], &
    [(bfloat16, integer :: j=1, count([bfloat16 >= 0]))], &
    [(real32, integer :: j=1, count([real32 >= 0]))], &
    [(real64, integer :: j=1, count([real64 >= 0]))], &
    [(real80, integer :: j=1, count([real80 >= 0]))], &
    [(real64x2, integer :: j=1, count([real64x2 >= 0]))], &
    [(real128, integer :: j=1, count([real128 >= 0]))]]

  integer, parameter :: current_team = -1, initial_team = -2, parent_team = -3

  integer, parameter :: input_unit = 5, output_unit = 6
  integer, parameter :: error_unit = 0
  integer, parameter :: iostat_end = -1, iostat_eor = -2
  integer, parameter :: iostat_inquire_internal_unit = -1

  integer, parameter :: character_storage_size = 8
  integer, parameter :: file_storage_size = 8
  integer, parameter :: numeric_storage_size = 32

  integer, parameter :: stat_failed_image = FORTRAN_RUNTIME_STAT_FAILED_IMAGE
  integer, parameter :: stat_locked = FORTRAN_RUNTIME_STAT_LOCKED
  integer, parameter :: stat_locked_other_image = FORTRAN_RUNTIME_STAT_LOCKED_OTHER_IMAGE
  integer, parameter :: stat_stopped_image = FORTRAN_RUNTIME_STAT_STOPPED_IMAGE
  integer, parameter :: stat_unlocked = FORTRAN_RUNTIME_STAT_UNLOCKED
  integer, parameter :: stat_unlocked_failed_image = FORTRAN_RUNTIME_STAT_UNLOCKED_FAILED_IMAGE

 contains

  character(len=80) function compiler_options()
    compiler_options = 'COMPILER_OPTIONS() not yet implemented'
  end function compiler_options

  character(len=80) function compiler_version()
    compiler_version = 'f18 in development'
  end function compiler_version
end module iso_fortran_env


text:
program main
  use xx_kinds
  use xx_prospino_subroutine
  implicit none

  integer                              :: inlo,isq_ng_in,icoll_in,i_error_in,ipart1_in,ipart2_in,isquark1_in,isquark2_in
  real(kind=double)                    :: energy_in
  integer                              :: ifs
  logical                              :: lfinal
  character(len=2)                     :: final_state_in

!----------------------------------------------------------------------------
  inlo = 1          ! specify LO only[0] or complete NLO (slower)[1]        !
!                   ! results: LO     - leading order, degenerate squarks   !
!                   !          NLO    - NLO, degenerate squarks             !
!                   !          LO_ms  - leading order, free squark masses   !
!                   !          NLO_ms - NLO, free squark masses             !
!                   ! all numerical errors (hopefully) better than 1%       !
!                   ! follow Vergas iteration on screen to check            !
!----------------------------------------------------------------------------

!----------------------------------------------------------------------------
  isq_ng_in = 1     ! specify degenerate [0] or free [1] squark masses      !
                    ! [0] means Prospino2.0 with average squark masses      !
                    ! [0] invalidates isquark_in switch                     !
!----------------------------------------------------------------------------

!----------------------------------------------------------------------------
  icoll_in = 1      ! collider : tevatron[0], lhc[1]                        !
  energy_in = 7000  ! collider energy in GeV                                !
!----------------------------------------------------------------------------

!----------------------------------------------------------------------------
  i_error_in = 0    ! with central scale [0] or scale variation [1]         !
!----------------------------------------------------------------------------

!----------------------------------------------------------------------------
!  final_state_in = ng   ! neutralino/chargino + gluino                     !
!                   ns   ! neutralino/chargino + squark                     !
!                   nn   ! neutralino/chargino pair combinations            !
!                   ll   ! slepton pair combinations                        !
!                   sb   ! squark-antisquark                                !
!                   ss   ! squark-squark                                    !
!                   tb   ! stop-antistop                                    !
!                   bb   ! sbottom-antisbottom                              !
!                   gg   ! gluino pair                                      !
!                   sg   ! squark + gluino                                  !
!                   lq   ! leptoquark pairs (using stop1 mass)              !
!                   le   ! leptoquark plus lepton (using stop1 mass)        !
!                   hh   ! charged Higgs pairs (private code only!)         !
!                   ht   ! charged Higgs with top (private code only!)      !
!                   xx   ! stop-antistop, wrapped Propino1                  !
!                                                                           !
!  squark and antisquark added, but taking into account different sb or ss  !
!----------------------------------------------------------------------------
  final_state_in = 'nn'

!----------------------------------------------------------------------------
!  final_state_in = ng,ns,nn                                                !
!  ipart1_in   = 1,2,3,4  neutralinos                                       !
!                5,6      positive charge charginos                         !
!                7,8      negative charge charginos                         !
!  ipart2_in the same                                                       !
!      chargino+ and chargino- different processes                          !
!                                                                           !
!  final_state_in = ll                                                      !
!  ipart1_in   = 0        sel,sel + ser,ser                                 !
!                1        sel,sel                                           !
!                2        ser,ser                                           !
!                3        snel,snel                                         !
!                4        sel+,snl                                          !
!                5        sel-,snl                                          !
!                6        stau1,stau1                                       !
!                7        stau2,stau2                                       !
!                8        stau1,stau2                                       !
!                9        sntau,sntau                                       !
!               10        stau1+,sntau                                      !
!               11        stau1-,sntau                                      !
!               12        stau2+,sntau                                      !
!               13        stau2-,sntau                                      !
!               14        H+,H- in Drell-Yan channel                        !
!                                                                           !
!  final_state_in = tb and bb                                               !
!  ipart1_in   = 1        stop1/sbottom1 pairs                              !
!                2        stop2/sbottom2 pairs                              !
!                                                                           !
!  for LO with light-squark flavor in the final state                       !
!  isquark1_in     =  -5,-4,-3,-2,-1,+1,+2,+3,+4,+5                         !
!                    (bL cL sL dL uL uR dR sR cR bR) in CteQ ordering       !
!  isquark1_in     = 0 sum over light-flavor squarks throughout             !
!                      (the squark mass in the data files is then averaged) !
!                                                                           !
!  flavors in initial state: only light-flavor partons, no bottoms          !
!                            bottom partons only for Higgs channels         !
!                                                                           !
!  flavors in final state: light-flavor quarks summed over five flavors     !
!                                                                           !
!----------------------------------------------------------------------------
  ipart1_in = 1
  ipart2_in = 1
  isquark1_in = 0
  isquark2_in = 0

  call PROSPINO_OPEN_CLOSE(0)                                                            ! open all input/output files

  do i_error_in=0,1,1                                                                    ! loop for test purposes only
  do icoll_in=0,1,1                                                                      ! loop for test purposes only
  do ifs = 1,15,1                                                                        ! loop and if structure for test purposes only
  if (ifs==1) then
     final_state_in = 'ng'
  else if (ifs==2) then
     final_state_in = 'ns'
  else if (ifs==3) then
     final_state_in = 'nn'
  else if (ifs==4) then
     final_state_in = 'll'
  else if (ifs==5) then
     final_state_in = 'gg'
  else if (ifs==6) then
     final_state_in = 'sb'
  else if (ifs==7) then
     final_state_in = 'sg'
  else if (ifs==8) then
     final_state_in = 'ss'
  else if (ifs==9) then
     final_state_in = 'tb'
  else if (ifs==10) then
     final_state_in = 'bb'
  else if (ifs==11) then
     final_state_in = 'lq'
  else if (ifs==12) then
     final_state_in = 'le'
  else if (ifs==13) then
     final_state_in = 'hh'
  else if (ifs==14) then
     final_state_in = 'ht'
  else if (ifs==15) then
     final_state_in = 'xx'
  end if

  do ipart1_in = 1,14,1                                                                   ! loop for test purposes only
  do ipart2_in = ipart1_in,6,1                                                                   ! loop for test purposes only

  call PROSPINO_CHECK_FS(final_state_in,ipart1_in,ipart2_in,lfinal)                      ! check final state
  if (.not. lfinal ) cycle

  call PROSPINO(inlo,isq_ng_in,icoll_in,energy_in,i_error_in,final_state_in,ipart1_in,ipart2_in,isquark1_in,isquark2_in) ! actual prospino call

  end do                                                                                ! close both ipart loops
  end do
  end do                                                                                ! close the ifs dummy loop
  end do                                                                                ! close the icoll_in loop
  end do                                                                                ! close the i_error_in loop

!----------------------------------------------------------------------------
!  input file: prospino.in.leshouches                                       !
!              use block MASS for masses, plus low-energy mixing matrices   !
!                                                                           !
!  output file: prospino.dat   for compact output format                    !
!               prospino.dat2  for long output including subchannels        !
!               prospino.dat3  lo file for masses, flags, etc               !
!----------------------------------------------------------------------------
  call PROSPINO_OPEN_CLOSE(1)                                                            ! close all input/output files

end program main

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


text:
! calmat.f for wcalprem.f
!----------------------------------------------------------
subroutine calmatc(nlayer,vnp,vra,con,rpow,w1dn,w2dn,ra,m,work)
!----------------------------------------------------------
! Computing \int r^rpow con W_p^(w1dn) W_q^(w2dn) dr.
!----------------------------------------------------------
    implicit none
    integer:: maxrpow
    parameter ( maxrpow = 2 )
    integer,intent(in):: nlayer,vnp,rpow,w1dn,w2dn
    double precision:: vra(vnp),con(vnp),ra(*),m(*),work(*)
    integer:: i,j,k,l,nn,snp
    double precision:: a(2,2),b(2,2),c(5),rh
    ! parameter check
    if (maxrpow<rpow) stop 'Invalid arguments.(calmatc)'
    ! computing of the matrix elements
    snp = 1
    do i=1,nlayer
        rh = ra(i+1) - ra(i)
        if ( w1dn == 0 ) then
            a(2,1) = -1.d0 / rh
            a(1,1) = ra(i+1) / rh
            a(2,2) = 1.d0 / rh
            a(1,2) = -ra(i) / rh
        elseif ( w1dn == 1 ) then
                a(2,1) = 0.d0
                a(1,1) = -1.d0 / rh
                a(2,2) = 0.d0
                a(1,2) = 1.d0 / rh
        else
            stop 'Invalid arguments.(calmatc)'
        endif

        if ( w2dn == 0 ) then
            b(2,1) = -1.d0 / rh
            b(1,1) = ra(i+1) / rh
            b(2,2) = 1.d0 / rh
            b(1,2) = -ra(i) / rh
        elseif ( w2dn == 1 ) then
                b(2,1) = 0.d0
                b(1,1) = -1.d0 / rh
                b(2,2) = 0.d0
                b(1,2) = 1.d0 / rh
        else
            stop 'Invalid arguments.(calmatc)'
        endif

        do j=1,2
            do k=1,2
                c(1:5) = 0.d0

                call pmulti( 2,a(1,j),2,b(1,k),3,c )
                do l=3,1,-1
                    c(l+rpow) = c(l)
                    if ( rpow>0 ) c(l)=0.d0
                enddo
                nn = 4 * (i-1) + 2 * (j-1) + k
                call pinteg(snp,5,c,ra(i),ra(i+1),vnp,vra,con,work(nn))
            enddo
        enddo
    enddo
    if ( w1dn /= w2dn ) then
        m(1:nlayer*4-3:4)=2*work(1:nlayer*4-3:4)
        m(2:nlayer*4-2:4)=work(2:nlayer*4-2:4)+work(3:nlayer*4-1:4)
        m(3:nlayer*4-1:4)=work(2:nlayer*4-2:4)+work(3:nlayer*4-1:4)
        m(4:nlayer*4:4)=2*work(4:nlayer*4:4)
    else
        m(1:4*nlayer) = work(1:4*nlayer)
    endif
    return
end

!----------------------------------------------------------
subroutine caltl( nlayer,vnp,vra,rho,ra,tl)
!----------------------------------------------------------
! Computing of lumped mass matrix.
!----------------------------------------------------------
    implicit none
    integer:: nlayer,vnp
    double precision:: vra(vnp),rho(vnp),ra(*),tl(*)
    integer:: i,nn,snp
    double precision:: c(3),from,to

    snp = 1
    do i=1,nlayer
        c(1:2) = 0.d0
        c(3) = 1.d0

        from = ra(i)
        to = ( ra(i) + ra(i+1) ) / 2.d0
        nn = 4 * (i-1)
        call pinteg( snp,3,c,from,to,vnp,vra,rho,tl(nn+1) )

        tl(nn+2:nn+3) = 0.d0

        from = to
        to = ra(i+1)
        call pinteg( snp,3,c,from,to,vnp,vra,rho,tl(nn+4) )

    enddo
    return
end

!----------------------------------------------------------
subroutine calhl( nlayer,vnp,vra,mu,ra,hl)
!----------------------------------------------------------
! Computing of lumped rigidity matrix.
!----------------------------------------------------------
    implicit none
    integer:: nlayer,vnp
    double precision:: vra(vnp),mu(vnp),ra(*),hl(*)
    integer:: i,nn,snp
    double precision:: c(1),from,to

    snp = 1
    do  i=1,nlayer
        c(1) = 1.d0

        from = ra(i)
        to = ( ra(i) + ra(i+1) ) / 2.d0
        nn = 4 * (i-1)
        call pinteg( snp,1,c,from,to,vnp,vra,mu,hl(nn+1) )

        hl(nn+2:nn+3) = 0.d0

        from = to
        to = ra(i+1)
        call pinteg( snp,1,c,from,to,vnp,vra,mu,hl(nn+4) )

    enddo
    return
end

!----------------------------------------------------------
subroutine pmulti(n,a,m,b,l,c)
!----------------------------------------------------------
! Computing the (l-1) degrees polynomial c(n) which is the product of
! (n-1) degrees polynomial a(n) and (m-1) degrees polynomial b(n).
!----------------------------------------------------------
    integer:: n,m,l,i,j
    double precision:: a(n),b(m),c(l)

    if (n+m-1/=l) stop 'Invalid arguments.(pmulti)'

    c(1:l)=0.d0

    do i=1,n
        do j=1,m
            c(i+j-1) = c(i+j-1) + a(i)*b(j)
        enddo
    enddo

    return
end

!----------------------------------------------------------
subroutine pinteg(snp,n,p,from,to,vnp,vra,con,pint)
!----------------------------------------------------------
! Evaluating the integrated value pint from 'from' to 'to' of p(n)*con
! which is the product of (n-1) degrees polynomial 'p(n)' and 'con'.
!----------------------------------------------------------
    implicit none
    integer:: maxn
    parameter ( maxn = 5 )

    integer:: snp,n,vnp
    double precision:: from,to,p(n),vra(vnp),con(vnp),pint
    double precision:: x1,x2,q(2),pq(maxn+1),psint

    if ( maxn<n )stop 'Degrees of a polynomial is too large.(pinteg)'
    if ( snp>=vnp ) snp = 1

    pint = 0.d0
    x1 = from
    do
        if (  vra(snp)<=x1 .and. vra(snp+1)>x1  ) then
            x2 = dmin1( to, vra(snp+1) )
        else
            snp = snp + 1
            if ( snp==vnp ) snp = 1
            cycle
        endif
        ! evaluating the integrated value
        q(2) = ( con(snp+1)-con(snp) ) / ( vra(snp+1)-vra(snp) )
        q(1) = - q(2) * vra(snp) + con(snp)
        call pmulti( n,p,2,q,n+1,pq )
        call polint( n+1,pq,x1,x2,psint )
        pint = pint + psint

        if ( x2==to ) exit
        x1 = x2
    enddo
    return
end

!----------------------------------------------------------
subroutine polint( n,p,x1,x2,pint )
!----------------------------------------------------------
!c Evaluating the integrated value 'pint' from 'x1' to 'x2'
!c of (n-1) degrees polynomial 'p(n)'.
!----------------------------------------------------------
    implicit none
    integer,parameter:: maxn=6
    integer:: n
    double precision:: p(*),x1,x2,pint
    integer:: i,j
    double precision:: a(maxn),b(maxn),dx,xx

    if (maxn<n) stop 'Degrees of a polynomial is too large.(polint)'

    a(1) = 1.d0
    b(1) = 1.d0
    if ( n>=2 ) then
        do i=2,n
            a(i) = a(i-1) * x1
            b(i) = b(i-1) * x2
        enddo
    endif
    dx = x2 - x1

    pint = 0.d0
    do i=1,n
        xx = 0.d0
        do j=1,i
            xx = xx + a(j) * b(i-j+1) / dble(i)
        enddo
        pint = pint + p(i) * dx * xx
    enddo
    return
end

!----------------------------------------------------------
subroutine overlap( nlayer,a,a2 )
!----------------------------------------------------------
! Overlapping the coefficient matrix elements in the solid part.
!----------------------------------------------------------
    use parameters
    implicit none
    integer,intent(in):: nlayer
    complex(dp):: a(*),a2(2,*)
    integer:: i,j,k,mu,m,i1,i2,k1,k2,nsize

    nsize = nlayer+1
    mu = 1
    m = mu + 1
    do j=1,nsize
        i1 = max0(1,j-mu)
        i2 = j
        do i=i1,i2
            k = i - j + m
            if ( i==j ) then
                if ( i==1 ) then
                    a2(k,j) = a2(k,j) + a(1)
                else
                    if ( i==nsize ) then
                        a2(k,j) = a2(k,j) + a( 4*nlayer )
                    else
                        k1 = 4 * i - 4
                        k2 = k1 + 1
                        a2(k,j) = a2(k,j) + a(k1) + a(k2)
                    endif
                endif
            endif
            if (i+1==j) then
                k1 = 4 * i - 2
                a2(k,j) = a2(k,j) + a(k1)
            endif
        enddo
    enddo
    return
end

!----------------------------------------------------------
subroutine calg2( l,m,spo,r0,mt,mu0,coef,ga,a,ga2,dr,g2 )
!----------------------------------------------------------
    use parameters
    implicit none

    integer:: l,m
    double precision:: spo,r0,mt(3,3),mu0
    complex(dp) coef,ga(*),a(*),ga2(2,*),g2(*)
    integer:: i,itmp
    double precision:: b1,sgn,eps,ier
    complex(dp) dd,cg2(3),dr(3),z(3)

    data eps/ -1.d0 /

    ! computing of particular solution with free surface boundary conditions
    cg2=0
    dd = 0
    if ( m>=0 ) then
        sgn = 1.d0
    else
        sgn = - 1.d0
    endif
    if ( iabs(m)==1 ) then
        b1 = dsqrt( dble( 2*l+1 ) / ( 16.d0 * pi ) )
        !c
        dd = dcmplx( b1 ) * ( dcmplx( sgn * mt(1,3), mt(1,2) ) )/ ( dcmplx( r0 * r0 * mu0 ) * coef )
        itmp = 4
        cg2=0
        do  i=2,3
            cg2(i) = - dd * ( ga(itmp+1) + ga(itmp+2) )
            itmp = itmp + 2
        enddo
    else
        if ( iabs(m)==2 ) then
            b1 = dsqrt( dble(2*l+1)*dble(l-1)*dble(l+2)/(64.d0*pi) )
            cg2(2) = dcmplx( b1 / r0 ) * dcmplx( 2.d0*mt(2,3), sgn*( mt(2,2)-mt(3,3) ) )
        endif
    endif
    if ( m==-2 .or. m==-l ) then
        call dclisb0( ga2,3,1,2,cg2,eps,dr,z,ier)
    else
        call dcsbsub0( ga2,3,1,2,cg2,eps,dr,z,ier)
    endif
    cg2(3) = cg2(3) + dd
    ! computation of the excitation vector
    itmp = dint(spo)
    g2(itmp+1) = a(1) * cg2(1) + a(2) * cg2(3)
    g2(itmp+2) = a(3) * cg2(1) + a(4) * cg2(3)

    return
end

max_stars_repo_name: seismobassoon/dsmpy

max_issues_repo_name: seismobassoon/dsmpy

max_forks_repo_name: seismobassoon/dsmpy


text:
module rans_boundary_average_advection
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ONE, TWO, HALF
    use mod_spalart_allmaras,   only: SA_sigma, SA_c_n1
    use mod_rans_efficient

    use type_operator,          only: operator_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use DNAD_D
    implicit none

    private

    !> Implementation of the Euler boundary average flux
    !!
    !!  - At a boundary interface, the solution states Q- and Q+ exists on opposite
    !!    sides of the boundary. The average flux is computed as Favg = 1/2(F(Q-) + F(Q+))
    !!
    !!  @author Nathan A. Wukie
    !!  @date   1/17/2019
    !!
    !--------------------------------------------------------------------------------
    type, extends(operator_t), public :: rans_boundary_average_advection_t

    contains

        procedure   :: init
        procedure   :: compute

    end type rans_boundary_average_advection_t
    !********************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   1/17/2019
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(rans_boundary_average_advection_t),   intent(inout) :: self

        ! Set operator name
        call self%set_name('RANS Boundary Average Advection')

        ! Set operator type
        call self%set_operator_type('Boundary Advective Flux')

        ! Set operator equations
        call self%add_primary_field('Density'   )
        call self%add_primary_field('Momentum-1')
        call self%add_primary_field('Momentum-2')
        call self%add_primary_field('Momentum-3')
        call self%add_primary_field('Energy'    )

    end subroutine init
    !********************************************************************************

    !>  Boundary Flux routine for Euler
    !!
    !!  @author Nathan A. Wukie
    !!  @date   1/17/2019
    !!
    !!-------------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(rans_boundary_average_advection_t),   intent(inout)   :: self
        type(chidg_worker_t),                               intent(inout)   :: worker
        class(properties_t),                                intent(inout)   :: prop

        ! Storage at quadrature nodes
        type(AD_D), allocatable, dimension(:) ::                            &
            density_m, mom1_m, mom2_m, mom3_m, energy_m,                    &
            density_p, mom1_p, mom2_p, mom3_p, energy_p,                    &
            u_m, v_m, w_m, invdensity_m, p_m, temperature_m, H_m,           &
            u_p, v_p, w_p, invdensity_p, p_p, temperature_p, H_p,           &
            flux_1_m, flux_2_m, flux_3_m,                                   &
            flux_1_p, flux_2_p, flux_3_p, r

        type(AD_D), allocatable, dimension(:)   ::  &
            density_nutilde_m, density_nutilde_p

        ! Interpolate solution to quadrature nodes
        density_m = worker%get_field('Density', 'value', 'face interior')
        density_p = worker%get_field('Density', 'value', 'face exterior')

        mom1_m    = worker%get_field('Momentum-1', 'value', 'face interior')
        mom1_p    = worker%get_field('Momentum-1', 'value', 'face exterior')

        mom2_m    = worker%get_field('Momentum-2', 'value', 'face interior')
        mom2_p    = worker%get_field('Momentum-2', 'value', 'face exterior')

        mom3_m    = worker%get_field('Momentum-3', 'value', 'face interior')
        mom3_p    = worker%get_field('Momentum-3', 'value', 'face exterior')

        energy_m  = worker%get_field('Energy', 'value', 'face interior')
        energy_p  = worker%get_field('Energy', 'value', 'face exterior')

        select case (trim(turbulence_model))
            case('Spalart Allmaras')
                density_nutilde_m = worker%get_field('Density * NuTilde', 'value', 'face interior')
                density_nutilde_p = worker%get_field('Density * NuTilde', 'value', 'face exterior')
            case('none')
                density_nutilde_m = ZERO*density_m
                density_nutilde_p = ZERO*density_p
            case default
                call chidg_signal_one(FATAL,"RANS Efficient: invalid turbulence_model.",trim(turbulence_model))
        end select

        ! Account for cylindrical. Get tangential momentum from angular momentum.
        if (worker%coordinate_system() == 'Cylindrical') then
            r = worker%coordinate('1','face interior')
            mom2_m = mom2_m / r
            mom2_p = mom2_p / r
        end if

        !
        ! Compute velocities
        !
        invdensity_m = ONE/density_m
        invdensity_p = ONE/density_p

        u_m = mom1_m*invdensity_m
        v_m = mom2_m*invdensity_m
        w_m = mom3_m*invdensity_m

        u_p = mom1_p*invdensity_p
        v_p = mom2_p*invdensity_p
        w_p = mom3_p*invdensity_p

        !
        ! Compute pressure
        !
        call compute_pressure_temperature(density_m,mom1_m,mom2_m,mom3_m,energy_m,p_m,temperature_m)
        call compute_pressure_temperature(density_p,mom1_p,mom2_p,mom3_p,energy_p,p_p,temperature_p)

        !
        ! Compute boundary condition energy and enthalpy
        !
        H_m = (energy_m + p_m)*invdensity_m
        H_p = (energy_p + p_p)*invdensity_p

        !----------------------------------
        !            mass flux
        !----------------------------------
        flux_1_m = (density_m * u_m)
        flux_2_m = (density_m * v_m)
        flux_3_m = (density_m * w_m)

        flux_1_p = (density_p * u_p)
        flux_2_p = (density_p * v_p)
        flux_3_p = (density_p * w_p)

        call worker%integrate_boundary_average('Density','Advection',       &
                                                flux_1_m,flux_2_m,flux_3_m, &
                                                flux_1_p,flux_2_p,flux_3_p)

        !----------------------------------
        !         momentum-1 flux
        !----------------------------------
        flux_1_m = (mom1_m * u_m) + p_m
        flux_2_m = (mom1_m * v_m)
        flux_3_m = (mom1_m * w_m)

        flux_1_p = (mom1_p * u_p) + p_p
        flux_2_p = (mom1_p * v_p)
        flux_3_p = (mom1_p * w_p)

        call worker%integrate_boundary_average('Momentum-1','Advection',    &
                                                flux_1_m,flux_2_m,flux_3_m, &
                                                flux_1_p,flux_2_p,flux_3_p)

        !----------------------------------
        !         momentum-2 flux
        !----------------------------------
        flux_1_m = (mom2_m * u_m)
        flux_2_m = (mom2_m * v_m) + p_m
        flux_3_m = (mom2_m * w_m)

        flux_1_p = (mom2_p * u_p)
        flux_2_p = (mom2_p * v_p) + p_p
        flux_3_p = (mom2_p * w_p)

        ! Convert to tangential to angular momentum flux
        if (worker%coordinate_system() == 'Cylindrical') then
            flux_1_m = flux_1_m * r
            flux_2_m = flux_2_m * r
            flux_3_m = flux_3_m * r

            flux_1_p = flux_1_p * r
            flux_2_p = flux_2_p * r
            flux_3_p = flux_3_p * r
        end if

        call worker%integrate_boundary_average('Momentum-2','Advection',    &
                                                flux_1_m,flux_2_m,flux_3_m, &
                                                flux_1_p,flux_2_p,flux_3_p)

        !----------------------------------
        !         momentum-3 flux
        !----------------------------------
        flux_1_m = (mom3_m * u_m)
        flux_2_m = (mom3_m * v_m)
        flux_3_m = (mom3_m * w_m) + p_m

        flux_1_p = (mom3_p * u_p)
        flux_2_p = (mom3_p * v_p)
        flux_3_p = (mom3_p * w_p) + p_p

        call worker%integrate_boundary_average('Momentum-3','Advection',    &
                                                flux_1_m,flux_2_m,flux_3_m, &
                                                flux_1_p,flux_2_p,flux_3_p)

        !----------------------------------
        !           energy flux
        !----------------------------------
        flux_1_m = (density_m * H_m * u_m)
        flux_2_m = (density_m * H_m * v_m)
        flux_3_m = (density_m * H_m * w_m)

        flux_1_p = (density_p * H_p * u_p)
        flux_2_p = (density_p * H_p * v_p)
        flux_3_p = (density_p * H_p * w_p)

        call worker%integrate_boundary_average('Energy','Advection',        &
                                                flux_1_m,flux_2_m,flux_3_m, &
                                                flux_1_p,flux_2_p,flux_3_p)

        !-----------------------------------------
        !            TURBULENCE FLUX
        !-----------------------------------------
        select case (trim(turbulence_model))
            case('Spalart Allmaras')
                flux_1_m = density_nutilde_m * u_m
                flux_2_m = density_nutilde_m * v_m
                flux_3_m = density_nutilde_m * w_m

                flux_1_p = density_nutilde_p * u_p
                flux_2_p = density_nutilde_p * v_p
                flux_3_p = density_nutilde_p * w_p
                call worker%integrate_boundary_average('Density * NuTilde','Advection', &
                                                        flux_1_m, flux_2_m, flux_3_m,   &
                                                        flux_1_p, flux_2_p, flux_3_p)
            case('none')

            case default
                call chidg_signal_one(FATAL,"RANS Efficient: invalid turbulence_model.",trim(turbulence_model))
        end select

    end subroutine compute
    !*********************************************************************************************************

end module rans_boundary_average_advection

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
Subroutine nnkaon(irun, iseed, ictrl, i1, i2, iblock, srt, pcx, pcy, pcz, nchrg)
  Parameter (maxstr=150001, maxr=1)
  Parameter (aka=0.498)
  Common /aa/r(3, maxstr)
  Common /bb/p(3, maxstr)
  Common /cc/e(maxstr)
  Common /ee/id(maxstr), lb(maxstr)
  Common /bg/betax, betay, betaz, gamma
  Common /nn/nnn
  Common /run/num
  Common /pa/rpion(3, maxstr, maxr)
  Common /pb/ppion(3, maxstr, maxr)
  Common /pc/epion(maxstr, maxr)
  Common /pd/lpion(maxstr, maxr)
  Dimension px(4), py(4), pz(4)
  Common /dpert/dpertt(maxstr, maxr), dpertp(maxstr), dplast(maxstr), dpdcy(maxstr), dpdpi(maxstr, maxr), dpt(maxstr, maxr), dpp1(maxstr, maxr), dppion(maxstr, maxr)
  Save
  dm3 = 0.938
  dm4 = 0.938
  n = 0
  If (nchrg<=-1 .Or. nchrg>=3) dm3 = 1.232
  If (nchrg==-2 .Or. nchrg==4) dm4 = 1.232
  iblock = 0
  Call fstate(iseed, srt, dm3, dm4, px, py, pz, iflag)
  If (iflag<0) Then
    ictrl = -1
    n = n + 1
    Return
  End If
  iblock = 12
  pxrota = px(1)
  pyrota = py(1)
  pzrota = pz(1)
  Call rotate(pcx, pcy, pcz, pxrota, pyrota, pzrota)
  px(1) = pxrota
  py(1) = pyrota
  pz(1) = pzrota
  pxrota = px(2)
  pyrota = py(2)
  pzrota = pz(2)
  Call rotate(pcx, pcy, pcz, pxrota, pyrota, pzrota)
  px(2) = pxrota
  py(2) = pyrota
  pz(2) = pzrota
  pxrota = px(3)
  pyrota = py(3)
  pzrota = pz(3)
  Call rotate(pcx, pcy, pcz, pxrota, pyrota, pzrota)
  px(3) = pxrota
  py(3) = pyrota
  pz(3) = pzrota
  pxrota = px(4)
  pyrota = py(4)
  pzrota = pz(4)
  Call rotate(pcx, pcy, pcz, pxrota, pyrota, pzrota)
  px(4) = pxrota
  py(4) = pyrota
  pz(4) = pzrota
  nnn = nnn + 2
  lpion(nnn, irun) = 23
  If (nchrg==-1 .Or. nchrg==-2) Then
  End If
  epion(nnn, irun) = aka
  lpion(nnn-1, irun) = 21
  epion(nnn-1, irun) = aka
  e1cm = sqrt(dm3**2+px(1)**2+py(1)**2+pz(1)**2)
  p1beta = px(1)*betax + py(1)*betay + pz(1)*betaz
  transf = gamma*(gamma*p1beta/(gamma+1)+e1cm)
  pt1i1 = betax*transf + px(1)
  pt2i1 = betay*transf + py(1)
  pt3i1 = betaz*transf + pz(1)
  eti1 = dm3
  lb1 = 2
  If (nchrg>=-2 .And. nchrg<=1) lb1 = 2
  If (nchrg==-2 .Or. nchrg==-1) Then
    lb1 = 6
  End If
  If (nchrg==1 .Or. nchrg==2) lb1 = 1
  If (nchrg==3 .Or. nchrg==4) lb1 = 9
  e2cm = sqrt(dm4**2+px(2)**2+py(2)**2+pz(2)**2)
  p2beta = px(2)*betax + py(2)*betay + pz(2)*betaz
  transf = gamma*(gamma*p2beta/(gamma+1)+e2cm)
  pt1i2 = betax*transf + px(2)
  pt2i2 = betay*transf + py(2)
  pt3i2 = betaz*transf + pz(2)
  eti2 = dm4
  lb2 = 2
  If (nchrg>=-1 .Or. nchrg<=1) lb2 = 2
  If (nchrg==2 .Or. nchrg==3) lb2 = 1
  If (nchrg==4) lb2 = 9
  If (nchrg==-2) lb2 = 6
  p(1, i1) = pt1i1
  p(2, i1) = pt2i1
  p(3, i1) = pt3i1
  e(i1) = eti1
  lb(i1) = lb1
  p(1, i2) = pt1i2
  p(2, i2) = pt2i2
  p(3, i2) = pt3i2
  e(i2) = eti2
  lb(i2) = lb2
  epcmk = sqrt(epion(nnn-1,irun)**2+px(3)**2+py(3)**2+pz(3)**2)
  betak = px(3)*betax + py(3)*betay + pz(3)*betaz
  transf = gamma*(gamma*betak/(gamma+1.)+epcmk)
  ppion(1, nnn-1, irun) = betax*transf + px(3)
  ppion(2, nnn-1, irun) = betay*transf + py(3)
  ppion(3, nnn-1, irun) = betaz*transf + pz(3)
  rpion(1, nnn-1, irun) = r(1, i1)
  rpion(2, nnn-1, irun) = r(2, i1)
  rpion(3, nnn-1, irun) = r(3, i1)
  dppion(nnn-1, irun) = dpertp(i1)*dpertp(i2)
  epcmak = sqrt(epion(nnn,irun)**2+px(4)**2+py(4)**2+pz(4)**2)
  betaak = px(4)*betax + py(4)*betay + pz(4)*betaz
  transf = gamma*(gamma*betaak/(gamma+1.)+epcmak)
  ppion(1, nnn, irun) = betax*transf + px(4)
  ppion(2, nnn, irun) = betay*transf + py(4)
  ppion(3, nnn, irun) = betaz*transf + pz(4)
  rpion(1, nnn, irun) = r(1, i2)
  rpion(2, nnn, irun) = r(2, i2)
  rpion(3, nnn, irun) = r(3, i2)
  dppion(nnn, irun) = dpertp(i1)*dpertp(i2)
  Return
End Subroutine nnkaon

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


PROGRAM main

  IMPLICIT NONE

  ! Pick a "long enough" string length
  INTEGER, PARAMETER :: string_len = 50
  INTEGER :: i, j, offset = 2
  ! Array holding alphabet and shifted alphabet
  CHARACTER(LEN=1), DIMENSION(26) :: letters, shifted_letters
  CHARACTER(LEN=string_len) :: word, ciphered_word

  !letters in order
  letters = (/"A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", &
      "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"/)

  PRINT*, "Enter word to cipher, in all caps:"
  !word should be all upper case and at most 30 characters
  READ(*, '(A)') word
  !word = "BISCUIT"

  !Create array shifted around first. This would work manually:
  !shifted_letters(1:26-offset) = letters(1+offset:26)
  !shifted_letters(26-offset:26) = letters(1:1+offset)

  ! OR we could also look for/remember a function to shift an array around
  shifted_letters = CSHIFT(letters, offset)

  ciphered_word = ""
  DO i = 1, string_len
    DO j = 1, 26
      ! Fortran doesn't let you refer to a single character in a string, but you can
      ! have a length-1 substring. I.e. you need start and end indices in slice
      IF(letters(j) == word(i:i)) THEN
        ciphered_word(i:i) = shifted_letters(j)
        ! Finish the inner loop
        EXIT
      END IF
    END DO
  END DO

  ! Fortran strings are exactly the length given - the end is padded with spaces
  ! so for output we should TRIM this off
  PRINT*, TRIM(word), " becomes ", TRIM(ciphered_word)

END PROGRAM main


text:
include 'sort_data_module.f90'

program ch2501
    ! Sorting Reals and Integers
    use iso_fortran_env
    use sort_data_module
    implicit none

    integer, parameter :: n = 1000000
    real(real64), allocatable, dimension(:) :: x
    integer(int32), allocatable, dimension(:) :: y
    integer :: allocate_status = 0

    print *, ' Program starts'

    allocate(x(1:n),stat=allocate_status)
    if (allocate_status/=0) then
        print *, ' Allocate failed.'
        print *, ' Program terminates'
        stop 10
    end if
    print *, ' Real allocate complete'
    call random_number(x)
    print *, ' Real array initialised'
    call sort_data(x,n)
    print *, ' Real sort ended'
    print *, ' First 10 reals'
    write(*,100) x(1:10)

    allocate(y(1:n),stat=allocate_status)
    if (allocate_status/=0) then
        print *, ' Allocate failed.'
        print *, ' Program terminates'
        stop 20
    end if
    y = int(x*1000000)
    deallocate(x)
    print *, ' Integer array initialised'
    call sort_data(y,n)
    print *, ' Sort ended'
    print *, ' First 10 integers'
    write(*,110) y(1:10)
    deallocate(y)
    print *, ' Deallocate'
    print *, ' Program terminates'

100 format (5(2x,e14.6))
110 format (5(2x,i10))
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
      program hessrst
      implicit none
c
c  $Id$
c  This program helps create a fdrst file for numerical hessian restarts.
c  There are several lines below that need to be modified so the program
c  picks up the correct hessian and fd_ddipole files and so that
c  some other info is set correctly.  Look for the comments.
c  First set the rank to be 3 times the number of atoms.
c
      integer rank
      parameter (rank = 141)
c
c     Restart file to be constructed
      CHARACTER*255 FILEATR
c     Ascii hessian (need for restart)
      CHARACTER*255 FILEHESS
c     derivative dipole file (need for restart)
      CHARACTER*255 FILEDDIPOLE

      double precision hess(rank,rank)
      double precision ddipole(3*rank)
      integer iatom_start, ixyz_start
      integer i, j, iflag_grad0, dipole_there

      logical does_it_exist
c
c
c  The following file names need to be set by the user.
c  FILEATR is the .fdrst file to be created, FILEHESS is
c  the .hessian file to be read from, and FILEDDIPOLE is
c  the .fd_ddipole file to be read from.
c
      FILEATR     = "test.fdrst"
      FILEHESS    = "test.hessian"
      FILEDDIPOLE = "test.fd_ddipole"
c
c  The following variables need to be set by the user.
c  iatom_start is the atom to start with, ixyz_start is
c  1, 2, or 3 for x, y, or z.
c
      iatom_start = 45
      ixyz_start  = 3
c
c  The user should not have to modify anything below here.
c

c
c Open the files
c
      does_it_exist = .false.
      inquire(file=FILEATR,exist=does_it_exist)
      if (does_it_exist) then
        write(6,*) 'There is already a .fdrst file.'
        stop
      else
        open(unit=69,file=FILEATR,
     &      form='unformatted',
     &      access='sequential',
     &      status='new')
      endif
c
      does_it_exist = .false.
      inquire(file=FILEHESS,exist=does_it_exist)
      if (does_it_exist) then
        open (unit=46, file=FILEHESS, form='formatted',
     &      access='sequential', status='unknown')
      else
        write(6,*) 'The hessian file does not exist.'
        stop
      endif
c
      does_it_exist = .false.
      inquire(file=FILEDDIPOLE,exist=does_it_exist)
      if (does_it_exist) then
        open (unit=67, file=FILEDDIPOLE, form='formatted',
     &      access='sequential', status='unknown')
      else
        write(6,*) 'The fd_ddipole file does not exist.'
        stop
      endif
c
c  Read hessian information
c
      rewind (unit=46)
      do i = 1, rank
        do j = 1, i
          read(46,*,end=6,err=6) hess(i,j)
          hess(j,i) = hess(i,j)
        enddo
      enddo
      close (unit=46,status='keep')
      write(6,*)
     &    '  Nuclear Hessian retrieved from ASCII file.'
c
c  Read dipole derivative information
c
      rewind (unit=67)
      read(67,*,end=6,err=6)(ddipole(i),i=1,3*rank)
      close (unit=67,status='keep')
      write(6,*)
     &   '  Dipole derivative retrieved from ASCII file.'
c
c  Write information to the fdrst file
c
      iflag_grad0 = 0
      dipole_there = 1
      write(69)iatom_start,ixyz_start,rank,iflag_grad0,dipole_there
      write(69)hess
      write(69)ddipole
      close(unit=69,status='keep')
      goto 7
c
    6 continue
      write(6,*) 'There is a problem reading a file'
    7 continue
c
      end

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
Coded by Stanley Rabinowitz, 12 Vine Brook Road, Westford MA, 01886-4212.
      INTEGER VECT(3350),BUFFER(201)
      DATA VECT/3350*2/,MORE/0/
      DO 2 N = 1,201
        KARRAY = 0
        DO 3 L = 3350,1,-1
          NUM = 100000*VECT(L) + KARRAY*L
          KARRAY = NUM/(2*L - 1)
    3     VECT(L) = NUM - KARRAY*(2*L - 1)
        K = KARRAY/100000
        BUFFER(N) = MORE + K
    2   MORE = KARRAY - K*100000
      WRITE (*,100) BUFFER
  100 FORMAT (I2,"."/(1X,10I5.5))
      END

max_stars_repo_name: mbirabhadra/RosettaCodeData


!-----------------------------------------------------------------------
!     detect active bands
!-----------------------------------------------------------------------
subroutine active_band
  use wf_module, only: occ
  use cpmd_variables,only:mstocck,MBC,MB,MB_0,MB_1,MBZ_0,MBZ_1,MSP_0,MSP_1 &
       ,MB_0_CPMD,MB_1_CPMD,ir_band_cpmd,id_band_cpmd
  use parallel_module, only: myrank,myrank_b,np_band
  implicit none
  integer :: s,k,n

  allocate( mstocck(MBZ_0:MBZ_1,MSP_0:MSP_1) )
  mstocck=0

  do s=MSP_0,MSP_1
  do k=MBZ_0,MBZ_1
  do n=1,MB
     if ( occ(n,k,s)>1.d-10 ) mstocck(k,s)=mstocck(k,s)+1
  enddo
  enddo
  enddo

  MBC=maxval(mstocck)

  if ( myrank == 0 ) then
     do s=MSP_0,MSP_1
     do k=MBZ_0,MBZ_1
        write(*,*) "occupation number of ",k,"-th BZ=",mstocck(k,s)
     enddo
     enddo
     write(*,*) "maximum occupation number of MB=",mbc
  endif

!--- for band-parallel calculation

  allocate( ir_band_cpmd(0:np_band-1) ) ; ir_band_cpmd= 0
  allocate( id_band_cpmd(0:np_band-1) ) ; id_band_cpmd=-1

  do n=1,MBC
     k=mod(n-1,np_band)
     ir_band_cpmd(k) = ir_band_cpmd(k) + 1
  end do
  do k=0,np_band-1
     id_band_cpmd(k) = sum( ir_band_cpmd(0:k) ) - ir_band_cpmd(k)
  end do

  MB_0_CPMD = id_band_cpmd(myrank_b) + 1
  MB_1_CPMD = id_band_cpmd(myrank_b) + ir_band_cpmd(myrank_b)

  return
end subroutine active_band


text:
module cuda_data_dispatcher_module
    use, intrinsic :: iso_c_binding, only : c_null_ptr, c_ptr

    use :: memory_api, only : &
            data_storage, &
            memory_manager
    use :: stream_module, only : stream
    use :: cuda_storage_module, only : cuda_storage

    implicit none
    private

    public :: cuda_data_dispatcher

    type :: cuda_data_dispatcher
    contains
        procedure :: get_temporary => get_temporary
        procedure :: update_from_temporary => update_from_temporary
        procedure :: release_temporary => release_temporary
        procedure :: update_and_release_temporary => update_and_release_temporary
    end type cuda_data_dispatcher

contains
    type(c_ptr) function get_temporary(this, storage, astream)
        class(cuda_data_dispatcher), intent(in) :: this
        class(data_storage), intent(in) :: storage
        type(stream), intent(in), optional :: astream

        type(cuda_storage) :: local

        local = cuda_storage()
        call memory_manager%secure_local_storage(storage, local, astream)
        get_temporary = local%get_data_pointer()

        call local%clear()
    end function get_temporary

    subroutine update_and_release_temporary(this, storage, src, astream)
        class(cuda_data_dispatcher), intent(in) :: this
        class(data_storage), intent(inout) :: storage
        type(c_ptr), intent(inout) :: src
        type(stream), intent(in), optional :: astream

        call this%update_from_temporary(storage, src, astream)
        call this%release_temporary(storage, src, astream)
    end subroutine update_and_release_temporary

    subroutine update_from_temporary(this, storage, src, astream)
        class(cuda_data_dispatcher), intent(in) :: this
        class(data_storage), intent(inout) :: storage
        type(c_ptr), intent(in) :: src
        type(stream), intent(in), optional :: astream

        type(cuda_storage) :: local_storage

        local_storage = cuda_storage(src, storage%number_of_bytes)
        call memory_manager%update_remote_storage(local_storage, storage, astream)
        call local_storage%clear()
    end subroutine update_from_temporary

    subroutine release_temporary(this, storage, ptr, astream)
        class(cuda_data_dispatcher), intent(in) :: this
        class(data_storage), intent(in) :: storage
        type(c_ptr), intent(inout) :: ptr
        type(stream), intent(in), optional :: astream

        type(cuda_storage) :: local_storage

        local_storage = cuda_storage(ptr, storage%number_of_bytes)
        call memory_manager%release_local_storage(storage, local_storage, astream)
        ptr = c_null_ptr
    end subroutine release_temporary
end module cuda_data_dispatcher_module

max_stars_repo_name: cheshyre/ntcl-data

max_issues_repo_name: cheshyre/ntcl-data

max_forks_repo_name: cheshyre/ntcl-data