!*==zgees.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief  ZGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download ZGEES + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
!                         LDVS, WORK, LWORK, RWORK, BWORK, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          JOBVS, SORT
!       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
!       ..
!       .. Array Arguments ..
!       LOGICAL            BWORK( * )
!       DOUBLE PRECISION   RWORK( * )
!       COMPLEX*16         A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
!       ..
!       .. Function Arguments ..
!       LOGICAL            SELECT
!       EXTERNAL           SELECT
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> ZGEES computes for an N-by-N complex nonsymmetric matrix A, the
!> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
!> vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
!>
!> Optionally, it also orders the eigenvalues on the diagonal of the
!> Schur form so that selected eigenvalues are at the top left.
!> The leading columns of Z then form an orthonormal basis for the
!> invariant subspace corresponding to the selected eigenvalues.
!>
!> A complex matrix is in Schur form if it is upper triangular.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] JOBVS
!> \verbatim
!>          JOBVS is CHARACTER*1
!>          = 'N': Schur vectors are not computed;
!>          = 'V': Schur vectors are computed.
!> \endverbatim
!>
!> \param[in] SORT
!> \verbatim
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the Schur form.
!>          = 'N': Eigenvalues are not ordered:
!>          = 'S': Eigenvalues are ordered (see SELECT).
!> \endverbatim
!>
!> \param[in] SELECT
!> \verbatim
!>          SELECT is a LOGICAL FUNCTION of one COMPLEX*16 argument
!>          SELECT must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'S', SELECT is used to select eigenvalues to order
!>          to the top left of the Schur form.
!>          IF SORT = 'N', SELECT is not referenced.
!>          The eigenvalue W(j) is selected if SELECT(W(j)) is true.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix A. N >= 0.
!> \endverbatim
!>
!> \param[in,out] A
!> \verbatim
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the N-by-N matrix A.
!>          On exit, A has been overwritten by its Schur form T.
!> \endverbatim
!>
!> \param[in] LDA
!> \verbatim
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!> \endverbatim
!>
!> \param[out] SDIM
!> \verbatim
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues for which
!>                         SELECT is true.
!> \endverbatim
!>
!> \param[out] W
!> \verbatim
!>          W is COMPLEX*16 array, dimension (N)
!>          W contains the computed eigenvalues, in the same order that
!>          they appear on the diagonal of the output Schur form T.
!> \endverbatim
!>
!> \param[out] VS
!> \verbatim
!>          VS is COMPLEX*16 array, dimension (LDVS,N)
!>          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
!>          vectors.
!>          If JOBVS = 'N', VS is not referenced.
!> \endverbatim
!>
!> \param[in] LDVS
!> \verbatim
!>          LDVS is INTEGER
!>          The leading dimension of the array VS.  LDVS >= 1; if
!>          JOBVS = 'V', LDVS >= N.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!> \endverbatim
!>
!> \param[in] LWORK
!> \verbatim
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,2*N).
!>          For good performance, LWORK must generally be larger.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!> \endverbatim
!>
!> \param[out] RWORK
!> \verbatim
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!> \endverbatim
!>
!> \param[out] BWORK
!> \verbatim
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value.
!>          > 0: if INFO = i, and i is
!>               <= N:  the QR algorithm failed to compute all the
!>                      eigenvalues; elements 1:ILO-1 and i+1:N of W
!>                      contain those eigenvalues which have converged;
!>                      if JOBVS = 'V', VS contains the matrix which
!>                      reduces A to its partially converged Schur form.
!>               = N+1: the eigenvalues could not be reordered because
!>                      some eigenvalues were too close to separate (the
!>                      problem is very ill-conditioned);
!>               = N+2: after reordering, roundoff changed values of
!>                      some complex eigenvalues so that leading
!>                      eigenvalues in the Schur form no longer satisfy
!>                      SELECT = .TRUE..  This could also be caused by
!>                      underflow due to scaling.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup complex16GEeigen
!
!  =====================================================================
      SUBROUTINE ZGEES(Jobvs,Sort,SELECT,N,A,Lda,Sdim,W,Vs,Ldvs,Work,   &
     &                 Lwork,Rwork,Bwork,Info)
      IMPLICIT NONE
!*--ZGEES201
!
!  -- LAPACK driver routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Jobvs , Sort
      INTEGER Info , Lda , Ldvs , Lwork , N , Sdim
!     ..
!     .. Array Arguments ..
      LOGICAL Bwork(*)
      DOUBLE PRECISION Rwork(*)
      COMPLEX*16 A(Lda,*) , Vs(Ldvs,*) , W(*) , Work(*)
!     ..
!     .. Function Arguments ..
      LOGICAL SELECT
      EXTERNAL SELECT
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ZERO , ONE
      PARAMETER (ZERO=0.0D0,ONE=1.0D0)
!     ..
!     .. Local Scalars ..
      LOGICAL lquery , scalea , wantst , wantvs
      INTEGER hswork , i , ibal , icond , ierr , ieval , ihi , ilo ,    &
     &        itau , iwrk , maxwrk , minwrk
      DOUBLE PRECISION anrm , bignum , cscale , eps , s , sep , smlnum
!     ..
!     .. Local Arrays ..
      DOUBLE PRECISION dum(1)
!     ..
!     .. External Subroutines ..
      EXTERNAL DLABAD , XERBLA , ZCOPY , ZGEBAK , ZGEBAL , ZGEHRD ,     &
     &         ZHSEQR , ZLACPY , ZLASCL , ZTRSEN , ZUNGHR
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      INTEGER ILAENV
      DOUBLE PRECISION DLAMCH , ZLANGE
      EXTERNAL LSAME , ILAENV , DLAMCH , ZLANGE
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MAX , SQRT
!     ..
!     .. Executable Statements ..
!
!     Test the input arguments
!
      Info = 0
      lquery = (Lwork==-1)
      wantvs = LSAME(Jobvs,'V')
      wantst = LSAME(Sort,'S')
      IF ( (.NOT.wantvs) .AND. (.NOT.LSAME(Jobvs,'N')) ) THEN
         Info = -1
      ELSEIF ( (.NOT.wantst) .AND. (.NOT.LSAME(Sort,'N')) ) THEN
         Info = -2
      ELSEIF ( N<0 ) THEN
         Info = -4
      ELSEIF ( Lda<MAX(1,N) ) THEN
         Info = -6
      ELSEIF ( Ldvs<1 .OR. (wantvs .AND. Ldvs<N) ) THEN
         Info = -10
      ENDIF
!
!     Compute workspace
!      (Note: Comments in the code beginning "Workspace:" describe the
!       minimal amount of workspace needed at that point in the code,
!       as well as the preferred amount for good performance.
!       CWorkspace refers to complex workspace, and RWorkspace to real
!       workspace. NB refers to the optimal block size for the
!       immediately following subroutine, as returned by ILAENV.
!       HSWORK refers to the workspace preferred by ZHSEQR, as
!       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
!       the worst case.)
!
      IF ( Info==0 ) THEN
         IF ( N==0 ) THEN
            minwrk = 1
            maxwrk = 1
         ELSE
            maxwrk = N + N*ILAENV(1,'ZGEHRD',' ',N,1,N,0)
            minwrk = 2*N
!
            CALL ZHSEQR('S',Jobvs,N,1,N,A,Lda,W,Vs,Ldvs,Work,-1,ieval)
            hswork = Work(1)
!
            IF ( .NOT.wantvs ) THEN
               maxwrk = MAX(maxwrk,hswork)
            ELSE
               maxwrk = MAX(maxwrk,N+(N-1)                              &
     &                  *ILAENV(1,'ZUNGHR',' ',N,1,N,-1))
               maxwrk = MAX(maxwrk,hswork)
            ENDIF
         ENDIF
         Work(1) = maxwrk
!
         IF ( Lwork<minwrk .AND. .NOT.lquery ) Info = -12
      ENDIF
!
      IF ( Info/=0 ) THEN
         CALL XERBLA('ZGEES ',-Info)
         RETURN
      ELSEIF ( lquery ) THEN
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N==0 ) THEN
         Sdim = 0
         RETURN
      ENDIF
!
!     Get machine constants
!
      eps = DLAMCH('P')
      smlnum = DLAMCH('S')
      bignum = ONE/smlnum
      CALL DLABAD(smlnum,bignum)
      smlnum = SQRT(smlnum)/eps
      bignum = ONE/smlnum
!
!     Scale A if max element outside range [SMLNUM,BIGNUM]
!
      anrm = ZLANGE('M',N,N,A,Lda,dum)
      scalea = .FALSE.
      IF ( anrm>ZERO .AND. anrm<smlnum ) THEN
         scalea = .TRUE.
         cscale = smlnum
      ELSEIF ( anrm>bignum ) THEN
         scalea = .TRUE.
         cscale = bignum
      ENDIF
      IF ( scalea ) CALL ZLASCL('G',0,0,anrm,cscale,N,N,A,Lda,ierr)
!
!     Permute the matrix to make it more nearly triangular
!     (CWorkspace: none)
!     (RWorkspace: need N)
!
      ibal = 1
      CALL ZGEBAL('P',N,A,Lda,ilo,ihi,Rwork(ibal),ierr)
!
!     Reduce to upper Hessenberg form
!     (CWorkspace: need 2*N, prefer N+N*NB)
!     (RWorkspace: none)
!
      itau = 1
      iwrk = N + itau
      CALL ZGEHRD(N,ilo,ihi,A,Lda,Work(itau),Work(iwrk),Lwork-iwrk+1,   &
     &            ierr)
!
      IF ( wantvs ) THEN
!
!        Copy Householder vectors to VS
!
         CALL ZLACPY('L',N,N,A,Lda,Vs,Ldvs)
!
!        Generate unitary matrix in VS
!        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
!        (RWorkspace: none)
!
         CALL ZUNGHR(N,ilo,ihi,Vs,Ldvs,Work(itau),Work(iwrk),           &
     &               Lwork-iwrk+1,ierr)
      ENDIF
!
      Sdim = 0
!
!     Perform QR iteration, accumulating Schur vectors in VS if desired
!     (CWorkspace: need 1, prefer HSWORK (see comments) )
!     (RWorkspace: none)
!
      iwrk = itau
      CALL ZHSEQR('S',Jobvs,N,ilo,ihi,A,Lda,W,Vs,Ldvs,Work(iwrk),       &
     &            Lwork-iwrk+1,ieval)
      IF ( ieval>0 ) Info = ieval
!
!     Sort eigenvalues if desired
!
      IF ( wantst .AND. Info==0 ) THEN
         IF ( scalea ) CALL ZLASCL('G',0,0,cscale,anrm,N,1,W,N,ierr)
         DO i = 1 , N
            Bwork(i) = SELECT(W(i))
         ENDDO
!
!        Reorder eigenvalues and transform Schur vectors
!        (CWorkspace: none)
!        (RWorkspace: none)
!
         CALL ZTRSEN('N',Jobvs,Bwork,N,A,Lda,Vs,Ldvs,W,Sdim,s,sep,      &
     &               Work(iwrk),Lwork-iwrk+1,icond)
      ENDIF
!
!
!        Undo balancing
!        (CWorkspace: none)
!        (RWorkspace: need N)
!
      IF ( wantvs ) CALL ZGEBAK('P','R',N,ilo,ihi,Rwork(ibal),N,Vs,Ldvs,&
     &                          ierr)
!
      IF ( scalea ) THEN
!
!        Undo scaling for the Schur form of A
!
         CALL ZLASCL('U',0,0,cscale,anrm,N,N,A,Lda,ierr)
         CALL ZCOPY(N,A,Lda+1,W,1)
      ENDIF
!
      Work(1) = maxwrk
!
!     End of ZGEES
!
      END SUBROUTINE ZGEES


text:
C
C $Id: draws.f,v 1.6 2008-07-27 00:17:25 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE DRAWS (MX1,MY1,MX2,MY2,IDRAW,IMARK)
C
C THIS ROUTINE DRAWS THE VISIBLE PART OF THE LINE CONNECTING
C (MX1,MY1) AND (MX2,MY2).  IF IDRAW .NE. 0, THE LINE IS DRAWN.
C IF IMARK .NE. 0, THE VISIBILITY ARRAY IS MARKED.
C
      LOGICAL         VIS1       ,VIS2
      DIMENSION       PXS(2)     ,PYS(2)
      COMMON /SRFBLK/ LIMU(1024) ,LIML(1024) ,CL(41)     ,NCL        ,
     1                LL         ,FACT       ,IROT       ,NDRZ       ,
     2                NUPPER     ,NRSWT      ,BIGD       ,UMIN       ,
     3                UMAX       ,VMIN       ,VMAX       ,RZERO      ,
     4                IOFFP      ,NSPVAL     ,SPVAL      ,BIGEST
      DATA STEEP/5./
      DATA MX,MY/0,0/
C
C  Do a call forcing a BLOCKDATA to be loaded from a binary library.
C
      CALL SRFABD
C
C MAKE LINE LEFT TO RIGHT.
C
      MMX1 = MX1
      MMY1 = MY1
      MMX2 = MX2
      MMY2 = MY2
      IF (MMX1.EQ.NSPVAL .OR. MMX2.EQ.NSPVAL) RETURN
      IF (MMX1 .GT. MMX2) GO TO  10
      NX1 = MMX1
      NY1 = MMY1
      NX2 = MMX2
      NY2 = MMY2
      GO TO  20
   10 NX1 = MMX2
      NY1 = MMY2
      NX2 = MMX1
      NY2 = MMY1
   20 IF (NUPPER .LT. 0) GO TO 180
C
C CHECK UPPER VISIBILITY.
C
      VIS1 = NY1 .GE. (LIMU(NX1)-1)
      VIS2 = NY2 .GE. (LIMU(NX2)-1)
C
C VIS1 AND VIS2 TRUE MEANS VISIBLE.
C
      IF (VIS1 .AND. VIS2) GO TO 120
C
C VIS1 AND VIS2 FALSE MEANS INVISIBLE.
C
      IF (.NOT.(VIS1 .OR. VIS2)) GO TO 180
C
C FIND CHANGE POINT.
C
      IF (NX1 .EQ. NX2) GO TO 110
      DY = REAL(NY2-NY1)/REAL(NX2-NX1)
      NX1P1 = NX1+1
      FNY1 = NY1
      IF (VIS1) GO TO  60
      DO  30 K=NX1P1,NX2
         MX = K
         MY = FNY1+REAL(K-NX1)*DY
         IF (MY .GT. LIMU(K)) GO TO  40
   30 CONTINUE
   40 IF (ABS(DY) .GE. STEEP) GO TO  90
   50 NX1 = MX
      NY1 = MY
      GO TO 120
   60 DO  70 K=NX1P1,NX2
         MX = K
         MY = FNY1+REAL(K-NX1)*DY
         IF (MY .LE. LIMU(K)) GO TO  80
   70 CONTINUE
   80 IF (ABS(DY) .GE. STEEP) GO TO 100
   81 NX2 = MX-1
      NY2 = MY
      GO TO 120
   90 IF (LIMU(MX) .EQ. 0) GO TO  50
      NX1 = MX
      NY1 = LIMU(NX1)
      GO TO 120
  100 IF (LIMU(MX-1).EQ.0) GO TO 81
      NX2 = MX-1
      NY2 = LIMU(NX2)
      GO TO 120
  110 IF (VIS1) NY2 = MIN(LIMU(NX1),LIMU(NX2))
      IF (VIS2) NY1 = MIN(LIMU(NX1),LIMU(NX2))
  120 IF (IDRAW .EQ. 0) GO TO 150
C
C DRAW VISIBLE PART OF LINE.
C
      IF (IROT) 130,140,130
  130 CONTINUE
      PXS(1) = REAL(NY1)
      PXS(2) = REAL(NY2)
      PYS(1) = REAL(1024-NX1)
      PYS(2) = REAL(1024-NX2)
      CALL GPL (2,PXS,PYS)
      GO TO 150
  140 CONTINUE
      PXS(1) = REAL(NX1)
      PXS(2) = REAL(NX2)
      PYS(1) = REAL(NY1)
      PYS(2) = REAL(NY2)
      CALL GPL (2,PXS,PYS)
  150 IF (IMARK .EQ. 0) GO TO 180
      IF (NX1 .EQ. NX2) GO TO 170
      DY = REAL(NY2-NY1)/REAL(NX2-NX1)
      FNY1 = NY1
      DO 160 K=NX1,NX2
         LTEMP = FNY1+REAL(K-NX1)*DY
         IF (LTEMP .GT. LIMU(K)) LIMU(K) = LTEMP
  160 CONTINUE
      GO TO 180
  170 LTEMP = MAX(NY1,NY2)
      IF (LTEMP .GT. LIMU(NX1)) LIMU(NX1) = LTEMP
  180 IF (NUPPER) 190,190,370
C
C SAME IDEA AS ABOVE, BUT FOR LOWER SIDE.
C
  190 IF (MMX1 .GT. MMX2) GO TO 200
      NX1 = MMX1
      NY1 = MMY1
      NX2 = MMX2
      NY2 = MMY2
      GO TO 210
  200 NX1 = MMX2
      NY1 = MMY2
      NX2 = MMX1
      NY2 = MMY1
  210 VIS1 = NY1 .LE. (LIML(NX1)+1)
      VIS2 = NY2 .LE. (LIML(NX2)+1)
      IF (VIS1 .AND. VIS2) GO TO 310
      IF (.NOT.(VIS1 .OR. VIS2)) GO TO 370
      IF (NX1 .EQ. NX2) GO TO 300
      DY = REAL(NY2-NY1)/REAL(NX2-NX1)
      NX1P1 = NX1+1
      FNY1 = NY1
      IF (VIS1) GO TO 250
      DO 220 K=NX1P1,NX2
         MX = K
         MY = FNY1+REAL(K-NX1)*DY
         IF (MY .LT. LIML(K)) GO TO 230
  220 CONTINUE
  230 IF (ABS(DY) .GE. STEEP) GO TO 280
  240 NX1 = MX
      NY1 = MY
      GO TO 310
  250 DO 260 K=NX1P1,NX2
         MX = K
         MY = FNY1+REAL(K-NX1)*DY
         IF (MY .GE. LIML(K)) GO TO 270
  260 CONTINUE
  270 IF (ABS(DY) .GE. STEEP) GO TO 290
  271 NX2 = MX-1
      NY2 = MY
      GO TO 310
  280 IF (LIML(MX) .EQ. 1024) GO TO 240
      NX1 = MX
      NY1 = LIML(NX1)
      GO TO 310
  290 IF (LIML(MX-1).EQ.1024) GO TO 271
      NX2 = MX-1
      NY2 = LIML(NX2)
      GO TO 310
  300 IF (VIS1) NY2 = MAX(LIML(NX1),LIML(NX2))
      IF (VIS2) NY1 = MAX(LIML(NX1),LIML(NX2))
  310 IF (IDRAW .EQ. 0) GO TO 340
      IF (IROT) 320,330,320
  320 CONTINUE
      PXS(1) = REAL(NY1)
      PXS(2) = REAL(NY2)
      PYS(1) = REAL(1024-NX1)
      PYS(2) = REAL(1024-NX2)
      CALL GPL (2,PXS,PYS)
      GO TO 340
  330 CONTINUE
      PXS(1) = REAL(NX1)
      PXS(2) = REAL(NX2)
      PYS(1) = REAL(NY1)
      PYS(2) = REAL(NY2)
      CALL GPL (2,PXS,PYS)
  340 IF (IMARK .EQ. 0) GO TO 370
      IF (NX1 .EQ. NX2) GO TO 360
      DY = REAL(NY2-NY1)/REAL(NX2-NX1)
      FNY1 = NY1
      DO 350 K=NX1,NX2
         LTEMP = FNY1+REAL(K-NX1)*DY
         IF (LTEMP .LT. LIML(K)) LIML(K) = LTEMP
  350 CONTINUE
      RETURN
  360 LTEMP = MIN(NY1,NY2)
      IF (LTEMP .LT. LIML(NX1)) LIML(NX1) = LTEMP
  370 RETURN
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
c======================================================================
      function bilinearintegral(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,
     &                                              z11,z12,z21,z22)
c======================================================================

      implicit none

*     !i/o
      double precision xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,
     &                                                z11,z12,z21,z22
      double precision bilinearintegral

*     !local
      double precision xlow,ylow,xhi,yhi,area,sumxi,sumeta,a,b,c,d

c#######################################################################
c
c         bilinearintegral integrates the bilinear with values z##
c         over the rectangular region xim <= x <= xip, and
c                                     yjm <= y <= yjp

c                                          written by David L George
c                                          Vancouver, WA April 2010
c#######################################################################

*     !integrate the portion of the bilinear intersected with the
c     !rectangular cell analytically.

c     !find limits of integral (this should already be true?)
      xlow = max(xim,x1)
      ylow = max(yjm,y1)
      xhi  = min(xip,x2)
      yhi =  min(yjp,y2)

*     !find the area of integration
      area = (yhi-ylow)*(xhi-xlow)
      sumxi = (xhi + xlow - 2.d0*x1)/dxx
      sumeta = (yhi + ylow - 2.d0*y1)/dyy

*     !find coefficients of bilinear a*xi + b*eta + c*xi*eta + d
      a = z21-z11
      b = z12-z11
      c = z22-z21-z12+z11
      d = z11

      bilinearintegral = (0.5d0*(a*sumxi + b*sumeta)
     &                     + 0.25d0*c*sumxi*sumeta + d)*area

      return

      end function

c======================================================================
      function bilinearintegral_s(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,
     &                                      z11,z12,z21,z22)
c======================================================================

      use geoclaw_module, only: r2d => rad2deg, d2r => deg2rad
      use geoclaw_module, only: Rearth => earth_radius

      implicit none

*     !i/o
      double precision xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,
     &                                 z11,z12,z21,z22
      double precision bilinearintegral_s

*     !local
      double precision xlo,ylo,xhi,yhi,delx,dely,a,b,c,d
      double precision xdiffhi,xdifflo,ydiffhi,ydifflo,xdiff2
      double precision adsinint,cbsinint

c#######################################################################
c
c         bilinearintegral integrates the bilinear with values z##
c         over the rectangular region xim <= x <= xip, and
c                                     yjm <= y <= yjp
c         integration is actually done on the surface of a sphere

c                                          written by David L George
c                                          Vancouver, WA April 2010
c#######################################################################

*     !integrate the portion of the bilinear intersected with the
c     !rectangular cell analytically.
c     !find limits of integral (this should already be true?)
      xlo = max(xim,x1)
      ylo = max(yjm,y1)
      xhi  = min(xip,x2)
      yhi =  min(yjp,y2)
      delx = xhi - xlo
      dely = yhi - ylo

*     !find terms for the integration
      xdiffhi = xhi - x1
      xdifflo = xlo - x1
      ydiffhi = yhi - y1
      ydifflo = ylo - y1
      xdiff2 = 0.5d0*(xdiffhi**2 - xdifflo**2)

      cbsinint = (r2d*cos(d2r*yhi) + ydiffhi*sin(d2r*yhi))
     &            -(r2d*cos(d2r*ylo) + ydifflo*sin(d2r*ylo))

      adsinint = r2d*(sin(d2r*yhi) - sin(d2r*ylo))

*     !find coefficients of bilinear a*xi + b*eta + c*xi*eta + d
      a = (z21-z11)/dxx
      b = (z12-z11)/dyy
      c = (z22-z21-z12+z11)/(dxx*dyy)
      d = z11

      bilinearintegral_s = ((a*xdiff2 + d*delx)*adsinint
     &            + r2d*(c*xdiff2 + b*delx)*cbsinint)*(Rearth*d2r)**2

      return

      end function

max_stars_repo_name: AsianHam/geoclaw

max_issues_repo_name: AsianHam/geoclaw

max_forks_repo_name: BrisaDavis/geoclaw


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2020 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine calcmach(vold,vcon,v,nk,
     &  ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,iout,
     &  nmethod,convergence,physcon,iponoel,inoel,ithermal,
     &  nactdoh,iit,compressible,ismooth,vcontu,vtu,turbulent,
     &  inomat,nodeboun,ndirboun,nboun,mi,co,factor)
!
!     calculates
!       vold (temperature,velocity and pressure)
!       vcon (volumetric energy density, volumetric momentum
!                density and density)
!       at the nodes
!
!     prints if iout=1
!
      implicit none
!
      integer convergence,compressible,
     &  nrhcon(*),ntmat_,nactdoh(0:4,*),iit,turbulent,mi(*),
     &  nshcon(*),ielmat(mi(3),*),nk,ithermal(*),i,j,k,index,iout,
     &  nmethod,imat,nelem,iponoel(*),inoel(3,*),ismooth,
     &  inomat(*),node,nodeboun(*),ndirboun(*),nboun
!
      real*8 v(0:mi(2),*),vold(0:mi(2),*),vcon(0:4,*),
     &  rhcon(0:1,ntmat_,*),rho,c1,vmax(0:4),dummy,press,
     &  voldmax(0:4),cp,r,temp,temp0,c2,c3,tempnew,vel2,
     &  shcon(0:3,ntmat_,*),drho,dtemp,physcon(*),dpress,
     &  vcontu(2,*),vtu(2,*),co(3,*),factor
!
!     calculate kappa (cp/cv) and store it in v(0,*)
!     calculate the Mach number and store it in v(1,*)
!
      do i=1,nk
         imat=inomat(i)
         temp=vold(0,i)
         call materialdata_cp_sec(imat,ntmat_,temp,shcon,
     &        nshcon,cp,physcon)
         r=shcon(3,1,imat)
         vel2=vold(1,i)**2+vold(2,i)**2+vold(3,i)**2
         v(0,i)=cp/(cp-r)
         v(1,i)=dsqrt((vold(1,i)**2+vold(2,i)**2+vold(3,i)**2)
     &        /(v(0,i)*r*(temp-physcon(1))))
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
Describe Users/BBrazill here.

Welcome to the Wiki! Thanks for getting involved with your El Macero Pharmacy business page! You might find Welcome to the Wiki/Business Owner this page helpful. Its not necessary to create a talk page just to add basic information to a page, so if theres any more info you want to add, go for it!  Users/ElleWeber

20091118 13:27:25 nbsp Wow Thank You  BBrazill      A young Man that works for you came in the other day. He had told me you loved them  Im glad  they work for you                                                        Thank You
                                                                     Steve
Users/LeeAnn

20110526 12:57:41 nbsp http://daviswiki.org/El_Macero_Pharmacy?actiondiff&version241&version140 Great response here!  Theres some work going on to rewrite the Welcome to the Wiki/Business Owner page right now, including examples of how to respond to complaints in a constructive and friendly manner.  Thanks providing a fantastic example! Users/TomGarberson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      PROGRAM TSCURVE
      IMPLICIT NONE
C
C  hycom_ts - Usage:  hycom_ts density
C
C                 generates a text file table of T and S at density
C                 for HYCOM's equation of state (sigma0 version).
C
C  this version for "serial" Unix systems.
C
C  Alan J. Wallcraft,  Naval Research Laboratory,  August 2001.
C
      INTEGER      IARGC
      INTEGER      NARG
      CHARACTER*240 CARG
C
      INTEGER      ITT
      REAL*8       RR,TT,SS
c
c --- coefficients for sigma-0 (based on Brydon & Sun fit)
      real*8     c1,c2,c3,c4,c5,c6,c7
      parameter (c1=-1.36471E-01, c2= 4.68181E-02, c3= 8.07004E-01,
     .           c4=-7.45353E-03, c5=-2.94418E-03,
     .           c6= 3.43570E-05, c7= 3.48658E-05)
c
c --- coefficients for sigma-2 (based on Brydon & Sun fit)
csig2 real       c1,c2,c3,c4,c5,c6,c7
csig2 parameter (c1= 9.77093E+00, c2=-2.26493E-02, c3= 7.89879E-01,
csig2.           c4=-6.43205E-03, c5=-2.62983E-03,
csig2.           c6= 2.75835E-05, c7= 3.15235E-05)
c
c --- salinity (mil) as a function of sigma and temperature (deg c)
      real*8 sofsig,r,t
      sofsig(r,t)=(r-c1-t*(c2+t*(c4+c6*t)))/(c3+t*(c5+c7*t))
C
C     READ ARGUMENTS.
C
      NARG = IARGC()
C
      IF     (NARG.EQ.1) THEN
        CALL GETARG(1,CARG)
        READ(CARG,*) RR
      ELSE
        WRITE(6,*) 'Usage: hycom_ts density'
        CALL EXIT(1)
      ENDIF
C
      WRITE(6,'(a,f7.3 / a)')
     &  '#   HYCOM equation of state, potential density =',RR,
     &  '#      T          S        RHO'
C
      DO ITT= -50,350
        TT = ITT*0.1
        SS = SOFSIG(RR,TT)
        WRITE(6,'(3f11.6)') TT,SS,RR
      ENDDO
      END

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
	SUBROUTINE SN_WPRT  ( isnfln, part, ihhmm, nlev, data, zwind,
     +			      iret )
C************************************************************************
C* SN_WPRT								*
C*									*
C* This subroutine writes data to an unmerged sounding data file.  The 	*
C* time and station must both be set before this subroutine is called.	*
C* The station time will be stored if the station time flag, STMFLG,	*
C* was set when the file was created.  This subroutine will only	*
C* write data to an unmerged data file.  The subroutine SN_WDAT must	*
C* be used to write data to a merged file.  The valid part names	*
C* are: TTAA, PPAA, TTBB, PPBB, TTCC, PPCC, TTDD, PPDD, TRPA, TRPC,     *
C* MXWA and MXWC.  The flag, ZWIND, is used only for significant wind   *
C* data (PPBB or PPDD).  If set, the winds are reported on height       *
C* surfaces; otherwise, the report is on pressure surfaces.	        *
C*									*
C* SN_WPRT  ( ISNFLN, PART, IHHMM, NLEV, DATA, ZWIND, IRET )		*
C*									*
C* Input parameters:							*
C*	ISNFLN		INTEGER		Sounding file number		*
C*	PART		CHAR*4		Part name			*
C*	IHHMM		INTEGER		Station time (HHMM)		*
C*	NLEV		INTEGER		Number of levels		*
C*	DATA (*)	REAL		Sounding data array		*
C*	ZWIND		LOGICAL		Flag for sig wind in z coord	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER	 	Return code			*
C*				    	   0 = normal return		*
C*				   	  -4 = file not open		*
C*					  -8 = station not set		*
C*				  	 -13 = DM error			*
C*					 -17 = invalid merge type	*
C*					 -22 = invalid part name	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 8/87						*
C* M. desJardins/GSFC	 9/87	Added significant winds on pres sfc.	*
C* M. desJardins/GSFC	10/87	Error in combining TTAA & PPAA data	*
C* M. desJardins/GSFC	 4/89	Added ABOVE in call to SN_CMAN		*
C* D. Kidwell/NCEP	 2/01	Removed mand. combine and merge, added  *
C*				check for tropopause part               *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
	INCLUDE		'sncmn.cmn'
C*
	CHARACTER*(*)	part
	REAL 		data (*)
	LOGICAL		zwind
C*
	CHARACTER	pprt*4
	LOGICAL		fix
	REAL		data2 ( 3, LLMXLV )
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	CALL SN_CHKF ( isnfln, iret )
	IF ( iret .ne. 0 )  RETURN
C
C*	Check that station has been set.
C
	irow = krow ( isnfln )
	icol = kcol ( isnfln )
	IF  ( ( irow .le. 0 ) .or. ( icol .le. 0 ) )  THEN
	    iret = -8
	    RETURN
	END IF
C
C*	Check that this is an unmerged data set.
C
	IF  ( mrgtyp ( isnfln ) )  THEN
	    iret = -17
	    RETURN
	END IF
C
C*	Compute number of data values to write.
C
	CALL ST_LCUC ( part, pprt, ier )
	IF  ( ( pprt .eq. 'TTAA' ) .or. ( pprt .eq. 'TTCC' ) )  THEN
	    ndata = nlev * 6
	  ELSE IF ( ( pprt .eq. 'TRPA' ) .or. ( pprt .eq. 'TRPC' )) THEN
	    ndata = nlev * 5
	  ELSE
	    ndata = nlev * 3
	END IF
C
C*	If this is significant wind data on a pressure surface, set the
C*	first valid pressure to be negative.
C
	fix = .false.
	IF  ( ( ( pprt .eq. 'PPBB' ) .or. ( pprt .eq. 'PPDD' ) ) .and.
     +	      ( .not. zwind ) )  THEN
	    fix = .true.
C
C*	    Move data to second array.
C
	    ipt = 1
	    DO  i = 1, nlev
		DO  j = 1, 3
		    data2 ( j, i ) = data ( ipt )
		    ipt = ipt + 1
		END DO
	    END DO
C
C*	    Make first valid level negative.
C
	    ilev = 1
	    DO WHILE  ( ilev .le. nlev )
		IF  ( .not. ERMISS ( data2 ( 1, ilev ) ) )  THEN
		    data2 ( 1, ilev ) = - data2 ( 1, ilev )
		    ilev = nlev + 1
		  ELSE
		    ilev = ilev + 1
		END IF
	    END DO
	END IF
C
C*	Write out the data.
C
	IF  ( fix )  THEN
	    CALL DM_WDTR  ( isnfln, irow, icol, pprt, ihhmm, data2,
     +			    ndata,  ier )
	  ELSE
	    CALL DM_WDTR  ( isnfln, irow, icol, pprt, ihhmm, data,
     +			    ndata,  ier )
	END IF
C
C*	Check for errors.
C
	IF  ( ier .eq. -10 )  THEN
	    iret = -22
	  ELSE IF  ( ier .ne. 0 )  THEN
	    iret = -13
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


      program demo_setprompts
      use M_kracken, only : kracken,iget,rget,sget,setprompts
      implicit none

      call setprompts('demo', ' -int My INTEGER value  -float My REAL value  -str My CHARACTER value')
      call kracken(   'demo', ' -int 100 -float 123.456 -str DEFAULT')
      write(*,'(a,i0)')'INTEGER IS ',iget('demo_int')
      write(*,'(a,g0)')'REAL IS ',rget('demo_float')
      write(*,'(a,a)')'STRING IS '//trim(sget('demo_str'))

      end program demo_setprompts


text:
 Bars
You can walk up to a bar just about anywhere and they can make you a screwdriver. If you want something a bit more unusual, its really the luck of the draw. Some drinks require ingredients that are not normally stocked. Other drinks have recipes that arent as widely known. List the more uncommon cocktails you can get in Davis here and tell us where you can find them.

  Our House offers specialty cocktails made from a number of lesscommon ingredients, such as houseinfused vodkas and freshlysqueezed juices. The menu changes seasonally and the one on the website is generally outdated.

  Tako Korean BBQ Tacos Tako is probably the only place in town to offer a cocktail made with Soju.

  At KetMoRee you can get a Full Moon Party Bowl, a concoction of rums and juices meant to be shared by 24 friends.

  GetDrunkQuick Schemes
Out for a night of adventure that youre not sure you want to remember? Several local bars have offerings to accommodate your needs.

  The Wicky Wacky Woo at Cafe Bernardo Bar Bernardo is the classic Davis drink for those who just need to get plastered.

  Red 88 Noodle Bar has its own take on the Wicky, called the FML.

  Many of the drinks on Sophias specialty cocktail menu are known to be rather strong, for example, the Thaibreaker.

 Where can I get.... ?
  New Orleans classics like the Ramos Gin Fizz and Sazerac?
    Sophias offers a Sazerac

 Recipes

If you want to make your own cocktails, the following pages may be of interest to you!

  Bloody Mary
  Mojito
  EctoKaze
  Eggnog

Is it still a cocktail if it is being served out of a water cooler? Perhaps, perhaps not. Regardless, these recipies are good for getting a lot of people drunk, fast.
  Jungle Juice
  Mizzle Fizzle Juice


!******************************************************************************
!
! Licensing:
!
!   This code is distributed under the MIT license.
!
! Modified:
!
!   22 April 2019
!
! Author:
!
!   Nescirem
!
!==============================================================================
!
!>M test_omphelloworld
!>  >S omphelloworld
!

    module test_omphelloworld

        use,intrinsic       :: iso_fortran_env, only: output_unit

        implicit none

        private
        public :: omphelloworld

    contains

        subroutine omphelloworld

        !openmp library
        use omp_lib

        implicit none

        integer             :: num_threads,tid,num_procs

        !parallel start and set 4 threads
        !$omp parallel default(private) num_threads(4)
        num_threads = omp_get_num_threads()
        tid = omp_get_thread_num()
        num_procs = omp_get_num_threads()

        write( output_unit,'(A,I2,A,I2,A,I2,A)' ) 'thread',tid,' of',num_threads,&
            ' threads is alive in',num_procs,' processes'
        !$omp end parallel

        end subroutine omphelloworld

    end module test_omphelloworld
!******************************************************************************

!******************************************************************************
#ifndef INTEGRATED_TESTS
    program test_omp_helloworld

        use test_omphelloworld , only: omphelloworld
        implicit none

        call omphelloworld

    end program test_omp_helloworld
#endif
!******************************************************************************


text:
!-------------------------------------------------------------------------
!Author: Matias Pablo, Borghi Orue
!Title: 2D Monte Carlo simulation of the Ising Model
!Version: 1
!Comments:  $ gfortran ranmar_seed.f95 subroutines.f95 kawasaki.f95 -o runv1
!INPUTS: None
!OUTPUTS: Two dat files: 'Magnetization30L.dat' and 'Energy30L.dat'
!-------------------------------------------------------------------------
program kawasaki

use subroutines
use rm
implicit none

	integer::nrows,ncols, average
	integer::posx,posy,dim
	integer:: MCS, EquilSteps, NumofTempPoints
	integer:: i,k,indexTemp !Dummy indeces
	integer :: index1, index2
	integer::exchangeSpin
	integer::Delta_E
	integer,allocatable::matriz(:,:)
	integer,allocatable::ip(:),im(:)
	integer::promediateCounter = 0
	real::rvec(1)
	real,allocatable::MagnetizationAver(:), MagnetizationAver2(:),MagnetizationAver4(:),EnergyAver(:),EnergyAver2(:)
	real,dimension(3)::DeltaE
	real::Tmax,Tmin,TempInterval,currentTemp,beta
	character(len=100)::MagFile ='Magnetization20L.dat', EnergFile = 'Energy20L.dat'

	! Read in input parameters from file "ising.in"
	open(unit=11,file="ConfigurationFile",status="old",action="read")
	read(11,*);read(11,*) nrows
	read(11,*);read(11,*) ncols
	read(11,*);read(11,*) MCS
	read(11,*);read(11,*) EquilSteps
	read(11,*);read(11,*) Tmax
	read(11,*);read(11,*) Tmin
	read(11,*);read(11,*) TempInterval
	read(11,*);read(11,*) average
	close(11)
	! Set the dimensions of the matrix of spin arrays. This program uses
	! periodic boundary conditions, so the first two rows and columns are
	! the same as the last two.
	allocate(matriz(nrows,ncols))

	allocate(MagnetizationAver(average))
	allocate(MagnetizationAver2(average))
	allocate(MagnetizationAver4(average))
	allocate(EnergyAver(average))
	allocate(EnergyAver2(average))

	!Change this to take in count that the matrix can be rectangular
	allocate(im(nrows))
	allocate(ip(nrows))

	do i = 1,nrows
		ip(i) = i+1
		im(i) = i-1
	end do
	ip(nrows) = 1
	im(1) = nrows

	dim = nrows*ncols
	MCS = dim*10
	EquilSteps = int (0.5*MCS)

	if ( EquilSteps >= MCS ) then
		STOP 'Change the value of EquilSteps in the Configuration File! It can not be greater than MCS'
	end if

	print *, "MCS = ",MCS, " EquilSteps = ", EquilSteps
	NumofTempPoints = int ( (Tmax - Tmin)/TempInterval ) + 1
	!If the file exists overwrite erase it.
	open (UNIT=1, FILE=MagFile,ACTION='write',STATUS='replace')
	write(1,*) "Temperature AverMagnetization AverMagnetization^2 Suceptibility Cumulant"
	open (UNIT=33, FILE=EnergFile,ACTION='write',STATUS='replace')
	write(33,*) "Temperature AverEnergy AverEnergy^2 C_v"

	!Initialize the matrix
	!with spin-up initial values
	call ini_matriz(matriz,nrows,ncols)

	do indexTemp = 1,NumofTempPoints !Loop over all the temperatures from Tmax to Tmin
		!Change the temperature value
		call New_Temp(indexTemp,currentTemp,Tmax,Tmin,TempInterval)
		beta = 1.0/currentTemp
		print *, "Current Temperature T = ", currentTemp

		!Initialize some variables to zero
		call SetToZeroVariables(MagnetizationAver,MagnetizationAver2, &
									MagnetizationAver4,EnergyAver,EnergyAver2,average)

		!Define the values of DeltaE knowing the Temperature
		!Only two possible values
		! 1) DeltaE(1) = exp(-4/Temp)
		! 2) DeltaE(2) = exp(-8/Temp)
		call Inicializacion_DeltaE(DeltaE,currentTemp)

		do k = 1 , average
			promediateCounter = 0
			!Choose a random point (posx,posy) where we begin calculating
			call elegir_sitio(nrows,ncols,posx,posy)

			!Sum over one MonteCarlo time //dim = 1MC
			!Sum over A MonteCarlo time: A*dim
			!In the next loop we flip a sping according to some conditions
			do i=1,MCS
				!Calculate the Magnetization and Energy after equilibrium
				if (i > EquilSteps ) then
					promediateCounter = promediateCounter + 1
					call CalculationsINequilibrium (nrows,ncols,matriz,Delta_E,EnergyAver,EnergyAver2, &
										MagnetizationAver,MagnetizationAver2,MagnetizationAver4,k,average,ip,im)

				end if
				!Calculate DeltaE for a given point (posx,posy)
				call CalcDeltaEnergy(nrows,ncols,posx,posy,matriz,Delta_E,ip,im)
				!print *, "Delta E: ",Delta_E
				! ..Flipping Spin..
				!If the new Energy spin-state is negative with flip the spin
				!If not it will depend on the probability exp(-DeltaE/Temp) > random point
				if (Delta_E <= 0) then
					exchangeSpin = matriz(posx,posy)
					if (posx == nrows) then
						matriz(posx,posy) = matriz(1,posy)
						matriz(1,posy) = exchangeSpin
					else
						matriz(posx,posy) = matriz(posx + 1,posy)
						matriz(posx+1,posy) = exchangeSpin
					end if
				else if(Delta_E == 4) then
					call ranmar(rvec,1)
					if( DeltaE(1) > rvec(1) ) then
						exchangeSpin = matriz(posx,posy)
						if (posx == nrows) then
							matriz(posx,posy) = matriz(1,posy)
							matriz(1,posy) = exchangeSpin
						else
							!print *,"OLD Pos x y: ", matriz(posx,posy)," Posx+1,y: ",matriz(posx+1,posy)
							matriz(posx,posy) = matriz(posx + 1,posy)
							matriz(posx+1,posy) = exchangeSpin
						end if
						!print *,"Spin exchanged for DeltaE = 1"
						!print *, "NEW Pos x y: ", matriz(posx,posy)," Posx+1,y: ",matriz(posx+1,posy)
					end if
				else if(Delta_E == 8) then
					call ranmar(rvec,1)
					if( DeltaE(2) > rvec(1) ) then
						exchangeSpin = matriz(posx,posy)
						if (posx == nrows) then
							matriz(posx,posy) = matriz(1,posy)
							matriz(1,posy) = exchangeSpin
						else
							matriz(posx,posy) = matriz(posx + 1,posy)
							matriz(posx+1,posy) = exchangeSpin
						end if
						!print *,"Spin exchanged for DeltaE = 6"
					end if
				else if(Delta_E == 12) then
					call ranmar(rvec,1)
					if( DeltaE(3) > rvec(1) ) then
						exchangeSpin = matriz(posx,posy)
						if (posx == nrows) then
							matriz(posx,posy) = matriz(1,posy)
							matriz(1,posy) = exchangeSpin
						else
							matriz(posx,posy) = matriz(posx + 1,posy)
							matriz(posx+1,posy) = exchangeSpin
						end if
						!print *,"Spin exchanged for DeltaE = 6"
					end if
				else
					STOP 'Error Delta_E take a value that is not one of the desired!!'
				end if

				!Choose a new random point (posx,posy) where we begin calculating
				call elegir_sitio(nrows,ncols,posx,posy)
				!----------------------------------------------------------------
				!-----PRINTF IN DYNAMIC LATTICE ---------------------------------
				!do index1=1,nrows
				!	do index2=1,ncols
				!		write(*,*) index1,index2,matriz(index1,index2)
				!	end do
				!end do
				!----------------------------------------------------------------
			end do !End loop of the MC steps
		end do ! End loop of averages
		!Write the results in a file
		call Write2File(promediateCounter, beta,currentTemp, MagnetizationAver,MagnetizationAver2, &
						MagnetizationAver4,EnergyAver,EnergyAver2,average)
	end do !End loop of Temperatures

	close(1)
	close(33)
end program kawasaki

max_stars_repo_path: V1 Kawasaki/kawasaki.f95

max_stars_repo_name: mattborghi/Thesis

max_issues_repo_path: V1 Kawasaki/kawasaki.f95

max_issues_repo_name: mattborghi/Thesis

max_forks_repo_path: V1 Kawasaki/kawasaki.f95

max_forks_repo_name: mattborghi/Thesis


MODULE mozart_mosaic_4bin_vbs0_Model

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Completely defines the model mozart_mosaic_4bin_vbs0
!    by using all the associated modules
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE mozart_mosaic_4bin_vbs0_Precision
  USE mozart_mosaic_4bin_vbs0_Parameters
  USE mozart_mosaic_4bin_vbs0_Global
  USE mozart_mosaic_4bin_vbs0_Function
  USE mozart_mosaic_4bin_vbs0_Integrator
  USE mozart_mosaic_4bin_vbs0_Rates
  USE mozart_mosaic_4bin_vbs0_Jacobian
  USE mozart_mosaic_4bin_vbs0_LinearAlgebra
  USE mozart_mosaic_4bin_vbs0_Monitor
  USE mozart_mosaic_4bin_vbs0_Util

END MODULE mozart_mosaic_4bin_vbs0_Model


text:
SUBROUTINE TEST09(m,n)
    REAL, DIMENSION(1:100,1:100) :: a,newa
    REAL w0,w1,w2,change,tolerance
    INTEGER i,j,iters

      change = tolerance + 1
      iters = 0
      do while ( change > tolerance )
         iters = iters + 1
         change = 0
!$xev dir append(copy(a),create(newa))
!$acc kernels
!$acc loop reduction(max:change)
         do j = 2, n-1
            do i = 2, m-1
               newa(i,j) = w0 * a(i,j) + &
               w1 * (a(i-1,j) + a(i,j-1) + a(i+1,j) + a(i,j+1) ) + &
               w2 * (a(i-1,j-1) + a(i-1,j+1) + a(i+1,j-1) + a(i+1,j+1) )
              change = max( change, abs( newa(i,j) - a(i,j) ) )
            end do
         end do
         a(2:m-1,2:n-1) = newa(2:m-1,2:n-1)
!$acc end kernels
      end do

    RETURN
END

max_stars_repo_name: xevolver/xevxml

max_issues_repo_name: xevolver/xevxml

max_forks_repo_name: xevolver/xevxml


! { dg-do compile }

program associate1
  type dl
    integer :: i
  end type
  type dt
    integer :: i
    real :: a(3, 3)
    type(dl) :: c(3, 3)
  end type
  integer :: v, i, j
  real :: a(3, 3)
  type(dt) :: b(3)
  i = 1
  j = 2
  associate(k => v, l => a(i, j), m => a(i, :))
  associate(n => b(j)%c(:, :)%i, o => a, p => b)
!$omp parallel shared (l)	! { dg-error "ASSOCIATE name" }
!$omp end parallel
!$omp parallel firstprivate (m)	! { dg-error "ASSOCIATE name" }
!$omp end parallel
!$omp parallel reduction (+: k)	! { dg-error "ASSOCIATE name" }
!$omp end parallel
!$omp parallel do firstprivate (k)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
!$omp parallel do lastprivate (n)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
!$omp parallel do private (o)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
!$omp parallel do shared (p)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
!$omp task private (k)		! { dg-error "ASSOCIATE name" }
!$omp end task
!$omp task shared (l)		! { dg-error "ASSOCIATE name" }
!$omp end task
!$omp task firstprivate (m)	! { dg-error "ASSOCIATE name" }
!$omp end task
!$omp do private (l)		! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
!$omp do reduction (*: k)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
!$omp sections private(o)	! { dg-error "ASSOCIATE name" }
!$omp section
!$omp section
!$omp end sections
!$omp parallel sections firstprivate(p)	! { dg-error "ASSOCIATE name" }
!$omp section
!$omp section
!$omp endparallelsections
!$omp parallelsections lastprivate(m)	! { dg-error "ASSOCIATE name" }
!$omp section
!$omp section
!$omp endparallelsections
!$omp sections reduction(+:k)	! { dg-error "ASSOCIATE name" }
!$omp section
!$omp section
!$omp end sections
!$omp simd private (l)		! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
  k = 1
!$omp simd lastprivate (m)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
  k = 1
!$omp simd reduction (+: k)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
  end do
  k = 1
!$omp simd linear (k : 2)	! { dg-error "ASSOCIATE name" }
  do i = 1, 10
    k = k + 2
  end do
  end associate
  end associate
end program


text:
      MODULE dot_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  08:33:41  03/09/06
      REAL(KIND(0.0D0)) FUNCTION dot (X, Y, N)
      USE vast_kind_param,ONLY: DOUBLE
      REAL(DOUBLE), DIMENSION(*), INTENT(IN) :: X
      REAL(DOUBLE), DIMENSION(*), INTENT(IN) :: Y
      INTEGER, INTENT(IN) :: N
      END FUNCTION
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
*DECK DAIGBT
      SUBROUTINE DAIGBT (RES, ADDA, NEQ, T, Y, YDOT,
     1                   MB, NB, PW, IPVT, IER )
      EXTERNAL RES, ADDA
      INTEGER NEQ, MB, NB, IPVT, IER
      INTEGER I, LENPW, LBLOX, LPB, LPC
      DOUBLE PRECISION T, Y, YDOT, PW
      DIMENSION Y(*), YDOT(*), PW(*), IPVT(*), NEQ(*)
C-----------------------------------------------------------------------
C This subroutine computes the initial value
C of the vector YDOT satisfying
C     A * YDOT = g(t,y)
C when A is nonsingular.  It is called by DLSOIBT for
C initialization only, when ISTATE = 0 .
C DAIGBT returns an error flag IER:
C   IER  =  0  means DAIGBT was successful.
C   IER .ge. 2 means RES returned an error flag IRES = IER.
C   IER .lt. 0 means the A matrix was found to have a singular
C              diagonal block (hence YDOT could not be solved for).
C-----------------------------------------------------------------------
      LBLOX = MB*MB*NB
      LPB = 1 + LBLOX
      LPC = LPB + LBLOX
      LENPW = 3*LBLOX
      DO 10 I = 1,LENPW
 10     PW(I) = 0.0D0
      IER = 1
      CALL RES (NEQ, T, Y, PW, YDOT, IER)
      IF (IER .GT. 1) RETURN
      CALL ADDA (NEQ, T, Y, MB, NB, PW(1), PW(LPB), PW(LPC) )
      CALL DDECBT (MB, NB, PW, PW(LPB), PW(LPC), IPVT, IER)
      IF (IER .EQ. 0) GO TO 20
      IER = -IER
      RETURN
 20   CALL DSOLBT (MB, NB, PW, PW(LPB), PW(LPC), YDOT, IPVT)
      RETURN
C----------------------- End of Subroutine DAIGBT ----------------------
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text: #redirect Davis Senior High School


! { dg-do compile }
program p
   integer :: n = 1
   character(3), parameter :: x(2) = ['abc', 'xyz']
   character(2), parameter :: y(2) = [x(2)(2:3), x(n)(1:2)] ! { dg-error "CHARACTER length must be a constant" }
end


text:
	SUBROUTINE IS_EGCA ( leng, report, iret )
C************************************************************************
C* IS_EGCA 								*
C*									*
C* This subroutine cleans up the point definition string of the area	*
C* covered by a MUHA international sigmet phenomenon:  Commas that	*
C* are used as decimal points in lat/lon values are replaced by		*
C* periods.  The remaining commas are replaced by blanks.  Slashes are  *
C* removed.								*
C*                                                                      *
C* IS_EGCA ( LENG, REPORT, IRET )					*
C*									*
C* Input parameters:							*
C*	LENG		INTEGER		Length of string                *
C*									*
C* Input and output parameters:						*
C*	REPORT		CHAR*		Partial sigmet report string    *
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* F. J. Yen/NCEP	10/01	Created					*
C************************************************************************
	CHARACTER*(*)	report
C------------------------------------------------------------------------
	iret = 0
C
C*	Replace commas used as decimal points with periods.
C*	Replace rest of commas with spaces.
C
	DO i = 1, leng - 3
	    IF ( report ( i:i ) .eq. ',' ) THEN
		IF ( report ( i + 2:i + 2 ) .eq. 'N' .or.
     +		     report ( i + 2:i + 2 ) .eq. 'W' .or.
     +		     report ( i + 2:i + 2 ) .eq. 'S' .or.
     +		     report ( i + 2:i + 2 ) .eq. 'E' )    THEN
		    CALL ST_ALNM ( report ( i - 1:i - 1 ), ityp, ier )
		    IF ( ityp .eq. 1) THEN
			report ( i:i ) = '.'
		      ELSE
			report ( i:i ) = ' '
		    END IF
		  ELSE IF ( report ( i + 3:i + 3 ) .eq. 'N' .or.
     +               	    report ( i + 3:i + 3 ) .eq. 'W' .or.
     +               	    report ( i + 3:i + 3 ) .eq. 'S' .or.
     +               	    report ( i + 3:i + 3 ) .eq. 'E' )    THEN
                    CALL ST_ALNM ( report ( i + 1:i + 1 ), ityp, ier )
                    IF ( ityp .eq. 1) THEN
                        CALL ST_ALNM ( report ( i + 2:i + 2 ),
     +				       ityp, ier )
                        IF ( ityp .eq. 1) THEN
			    report ( i:i ) = '.'
			  ELSE
			    report ( i:i ) = ' '
			END IF
                      ELSE
                        report ( i:i ) = ' '
		    END IF
		  ELSE
		    report ( i:i ) = ' '
		END IF
	    END IF
	END DO
C
C*	Remove slashes.
C
	ipos = - 1
	ib = 3
	ie = leng - 3
	DO WHILE ( ipos .ne. 0 .and. ie .gt. ib )
	    ipos = INDEX ( report ( ib:ie ), '/' )
	    IF ( ipos .ne. 0 ) THEN
		report = report ( 1:ib + ipos - 2 ) //
     +			    report ( ib + ipos : ie ) // ' '
		ie = ie - 1
		ib = ib + ipos - 1
	    END IF
	END DO
C*
	RETURN
	END


text:
      subroutine gf_unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,
     &                      mappdslen,coordlist,numcoord,ierr)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    gf_unpack4
!   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
!
! ABSTRACT: This subroutine unpacks Section 4 (Product Definition Section)
!   starting at octet 6 of that Section.
!
! PROGRAM HISTORY LOG:
! 2000-05-26  Gilbert
! 2002-01-24  Gilbert  - Changed to dynamically allocate arrays
!                        and to pass pointers to those arrays through
!                        the argument list.
!
! USAGE:    CALL gf_unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,mappdslen,
!    &                   coordlist,numcoord,ierr)
!   INPUT ARGUMENT LIST:
!     cgrib    - Character array that contains the GRIB2 message
!     lcgrib   - Length (in bytes) of GRIB message array cgrib.
!     iofst    - Bit offset of the beginning of Section 4.
!
!   OUTPUT ARGUMENT LIST:
!     iofst    - Bit offset of the end of Section 4, returned.
!     ipdsnum  - Product Definition Template Number ( see Code Table 4.0)
!     ipdstmpl - Pointer to integer array containing the data values for
!                the specified Product Definition
!                Template ( N=ipdsnum ).  Each element of this integer
!                array contains an entry (in the order specified) of Product
!                Defintion Template 4.N
!     mappdslen- Number of elements in ipdstmpl().  i.e. number of entries
!                in Product Defintion Template 4.N  ( N=ipdsnum ).
!     coordlist- Pointer to real array containing floating point values
!                intended to document
!                the vertical discretisation associated to model data
!                on hybrid coordinate vertical levels.  (part of Section 4)
!     numcoord - number of values in array coordlist.
!     ierr     - Error return code.
!                0 = no error
!                5 = "GRIB" message contains an undefined Product Definition
!                    Template.
!                6 = memory allocation error
!
! REMARKS: Uses Fortran 90 module pdstemplates and module re_alloc.
!
! ATTRIBUTES:
!   LANGUAGE: Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      use pdstemplates
      use re_alloc        !  needed for subroutine realloc

      character(len=1),intent(in) :: cgrib(lcgrib)
      integer,intent(in) :: lcgrib
      integer,intent(inout) :: iofst
      real,pointer,dimension(:) :: coordlist
      integer,pointer,dimension(:) :: ipdstmpl
      integer,intent(out) :: ipdsnum
      integer,intent(out) :: ierr,numcoord

      real(4),allocatable :: coordieee(:)
      integer,allocatable :: mappds(:)
      integer :: mappdslen
      logical needext

      ierr=0
      nullify(ipdstmpl,coordlist)

      call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
      iofst=iofst+32
      iofst=iofst+8     ! skip section number
      allocate(mappds(lensec))

      call gbyte(cgrib,numcoord,iofst,16)    ! Get num of coordinate values
      iofst=iofst+16
      call gbyte(cgrib,ipdsnum,iofst,16)    ! Get Prod. Def Template num.
      iofst=iofst+16
      !   Get Product Definition Template
      call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret)
      if (iret.ne.0) then
        ierr=5
        if( allocated(mappds) ) deallocate(mappds)
        return
      endif
      !
      !   Unpack each value into array ipdstmpl from the
      !   the appropriate number of octets, which are specified in
      !   corresponding entries in array mappds.
      !
      istat=0
      if (mappdslen.gt.0) allocate(ipdstmpl(mappdslen),stat=istat)
      if (istat.ne.0) then
         ierr=6
         nullify(ipdstmpl)
         if( allocated(mappds) ) deallocate(mappds)
         return
      endif
      do i=1,mappdslen
        nbits=iabs(mappds(i))*8
        if ( mappds(i).ge.0 ) then
          call gbyte(cgrib,ipdstmpl(i),iofst,nbits)
        else
          call gbyte(cgrib,isign,iofst,1)
          call gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
          if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
        endif
        iofst=iofst+nbits
      enddo
      !
      !   Check to see if the Product Definition Template needs to be
      !   extended.
      !   The number of values in a specific template may vary
      !   depending on data specified in the "static" part of the
      !   template.
      !
      if ( needext ) then
        call extpdstemplate(ipdsnum,ipdstmpl,newmappdslen,mappds)
        call realloc(ipdstmpl,mappdslen,newmappdslen,istat)
        !   Unpack the rest of the Product Definition Template
        do i=mappdslen+1,newmappdslen
          nbits=iabs(mappds(i))*8
          if ( mappds(i).ge.0 ) then
            call gbyte(cgrib,ipdstmpl(i),iofst,nbits)
          else
            call gbyte(cgrib,isign,iofst,1)
            call gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
            if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
          endif
          iofst=iofst+nbits
        enddo
        mappdslen=newmappdslen
      endif
      if( allocated(mappds) ) deallocate(mappds)
      !
      !   Get Optional list of vertical coordinate values
      !   after the Product Definition Template, if necessary.
      !
      nullify(coordlist)
      if ( numcoord .ne. 0 ) then
         allocate (coordieee(numcoord),stat=istat1)
         allocate(coordlist(numcoord),stat=istat)
         if ((istat1+istat).ne.0) then
            ierr=6
            nullify(coordlist)
            if( allocated(coordieee) ) deallocate(coordieee)
            return
         endif
        call gbytes(cgrib,coordieee,iofst,32,0,numcoord)
        call rdieee(coordieee,coordlist,numcoord)
        deallocate (coordieee)
        iofst=iofst+(32*numcoord)
      endif

      return    ! End of Section 4 processing
      end

max_stars_repo_name: martinremy/wps

max_issues_repo_name: martinremy/wps

max_forks_repo_name: WRF-CMake/WPS


text:
      subroutine rblankc(ialf,nch)   !  remove ALL blanks
      character ialf*(*),ich
      j=0
      nchold=nch
      do 10 i=1,nch
      ich=ialf(i:i)
      if(ich.eq.' ') go to 10
      j=j+1
      ialf(j:j)=ich
 10   continue
      nch=j
      if(nch.lt.nchold) ialf(nch+1:nchold)=' '
      return
      end

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
      PROGRAM xtutest
C     driver for routine tutest
      INTEGER NPTS,NSHFT,MPTS
      REAL EPS,VAR1,VAR2
      PARAMETER(NPTS=5000,MPTS=1000,EPS=0.02,VAR1=1.0,
     *     VAR2=4.0,NSHFT=10)
      INTEGER i,idum,j
      REAL data1(NPTS),data2(MPTS),fctr1,fctr2,gasdev,prob,shift,t
C     generate two Gaussian distributions of different variance
      idum=-51773
      fctr1=sqrt(VAR1)
      do 11 i=1,NPTS
        data1(i)=fctr1*gasdev(idum)
11    continue
      fctr2=sqrt(VAR2)
      do 12 i=1,MPTS
        data2(i)=(NSHFT/2.0)*EPS+fctr2*gasdev(idum)
12    continue
      write(*,'(1x,a,f6.2)') 'Distribution #1 : variance = ',VAR1
      write(*,'(1x,a,f6.2/)') 'Distribution #2 : variance = ',VAR2
      write(*,'(1x,t4,a,t18,a,t25,a)') 'Shift','T','Probability'
      do 14 i=1,NSHFT+1
        call tutest(data1,NPTS,data2,MPTS,t,prob)
        shift=(i-1)*EPS
        write(*,'(1x,f6.2,2f12.2)') shift,t,prob
        do 13 j=1,NPTS
          data1(j)=data1(j)+EPS
13      continue
14    continue
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xtutest.for


text:
Describe Users/austinchang here.

This persons real name is actually Yee Yang.


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE MCPIN (NAME,IK,NTESTG,INPCI)
!                                                                      *
!   This routine controls the computation  and storage of the values   *
!   and all indices of the angular coefficients                        *
!                                                                      *
!                                                                      *
!                   T  (ab)                                            *
!                    rs                                                *
!                                                                      *
!   k is the multipolarity of a two-particle Coulomb integral. a, b,   *
!   are orbital sequence numbers.  r and s are configuration           *
!   state function indices.                                            *
!                                                                      *
!   Call(s) to: [LIB92]: ALLOC, DALLOC, TNSRJJ                         *
!               [GENMCP]: QSORT.                                       *
!                                                                      *
!   Written by Per Jonsson                Last revision: JUne 1996     *
!   Bug corrected 2005-10-18                                           *
!                                                                      *
!***********************************************************************
!...Translated by Charlotte Froese Fischer
!                       Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE parameter_def,   ONLY: KEYORB, NNNW
      USE memory_man
      USE def_C
      USE eigv_C
      USE foparm_C
      USE mcp_C
      USE orb_C,           ONLY: ncf, nw
      USE prnt_C
      USE syma_C
      USE stat_C
      USE jqjc_C
      USE orbord_C
      USE mcpdata_C
      USE sbdat_C,         ONLY: KAMAX,NSHLII,NSHLFF,NSHLPPII,NSHLPPFF,&
                                 NINII,NINFF,NSHLPII,NSHLPFF,NAKINVII, &
                                 NAKINVFF,NLMAX
      USE sbdat1_C,        ONLY: nshlp, nshlpp
      USE cimat_C,         ONLY: cici, cfci
      USE blk_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
!      USE tnsrjj_I
!      USE cor_I
!      USE cord_I
!      USE itjpo_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      CHARACTER(LEN=24), INTENT(IN) :: name
      INTEGER, INTENT(IN) :: ik, ntestg, inpci
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!   Set the encoding key and its square
      INTEGER, PARAMETER :: KEY = KEYORB, KEYSQ = KEY*KEY, nvmax=100
      REAL(DOUBLE), PARAMETER :: CUTOFF = 1.0D-10
      INTEGER, PARAMETER :: NF = 200

      LOGICAL ::  F0INT,LINCR,RESTRT,COMP,AVAIL
      CHARACTER(LEN=20) :: CNUM
      CHARACTER(LEN=2)  :: CK
!
      REAL(DOUBLE), DIMENSION(NNNW) :: tshell
      REAL(DOUBLE), DIMENSION(NVMAX) ::SC
      INTEGER, DIMENSION(NNNW) :: nakinv
      INTEGER, DIMENSION(NLMAX) :: llistt, nshl, ninl
      REAL(DOUBLE), DIMENSION(20*NLMAX*NLMAX) :: cirot
      REAL(DOUBLE), DIMENSION(10000) :: EVSC
!
      REAL(DOUBLE), DIMENSION(:), pointer :: scr, ciout
      INTEGER :: ntest, ntestL, i, j, iblk, ncfd, nwd, kamaxd,  &
                 l, iioff, nvecsize
!
      NTESTL = 00
      NTEST = MAX(NTESTG,NTESTL)
      NTEST = 00
!
!   Set up data for the initial and final state case respectively
!
      IF (IK.EQ.1) THEN
        DO I = 1,NLMAX
          NSHL(I) = NSHLII(I)
          NINL(I) = NINII(I)
        ENDDO

        DO I = 1,NNNW
          NAKINV(I) = NAKINVII(I)
        ENDDO

        DO J = 1,NLMAX
          DO I = 1,NLMAX
            NSHLP(I,J) = NSHLPII(I,J)
          ENDDO
        ENDDO

        DO J = 1,NNNW
          DO I = 1,NLMAX
            NSHLPP(I,J) = NSHLPPII(I,J)
          ENDDO
        ENDDO

        DO I = 1,20*NLMAX*NLMAX
          CIROT(I) = CICI(I)
        ENDDO

      ELSEIF (IK.EQ.2) THEN

        DO I = 1,NLMAX
          NSHL(I) = NSHLFF(I)
          NINL(I) = NINFF(I)
        ENDDO

        DO I = 1,NNNW
          NAKINV(I) = NAKINVFF(I)
        ENDDO

        DO J = 1,NLMAX
          DO I = 1,NLMAX
            NSHLP(I,J) = NSHLPFF(I,J)
          ENDDO
        ENDDO

        DO J = 1,NNNW
          DO I = 1,NLMAX
            NSHLPP(I,J) = NSHLPPFF(I,J)
          ENDDO
        ENDDO

        DO I = 1,20*NLMAX*NLMAX
          CIROT(I) = CFCI(I)
        ENDDO
      ENDIF

      REWIND (NF)
      DO 1000 IBLK = 1, NBLOCK
!
!   Read the CI vectors
!
      CALL GETMIX(NAME,INPCI,IBLK)
!
!   Allocate memory for use in CITRAG
!
      CALL alloc (scr, ncf*nvec, 'scr', 'MCPIN')
      CALL alloc (ciout, ncf*nvec, 'CIOUT', 'MCPIN')
!
!   if data available on file read the datafile. Else sort the
!   MCP data into inegral based lists.
!
      READ(NF) NCFD,NWD,KAMAXD
      DO L = 1,KAMAX

!************
!
!. Offset for given L in shell matrices
!
          IF (L.EQ.1) THEN
            IIOFF = 1
          ELSE
            IIOFF = IIOFF + NSHL(L-1)** 2
          END IF
!  Corrected PER J

!**************

          READ(NF) NINTG,NCOEFF
          IF(NCOEFF*NINTG.EQ.0) CYCLE
!
!   Allocate memory. Note that in this case we can,
!   since the data is sorted, allocate less
!   memory than is done in qqsort.
!
          CALL ALLOC (JANN,NCOEFF, 'JANN', 'MCPIN')
          CALL ALLOC (JBNN,NCOEFF, 'JBNN', 'MCPIN')
          CALL ALLOC (INTGRL,NINTG, 'INTGRL', 'MCPIN')
          CALL ALLOC (CNN,NCOEFF, 'CNN', 'MCPIN')
          CALL ALLOC (INTPTR,NINTG, 'INTPTR', 'MCPIN')

          DO I = 1,NINTG
            READ(NF) INTGRL(I),INTPTR(I)
          ENDDO
          DO I = 1,NCOEFF
            READ(NF) CNN(I),JANN(I),JBNN(I)
          ENDDO
!*
!*. Offset for given L in shell matrices
!*
!        IF (L.EQ.1) THEN
!          IIOFF = 1
!        ELSE
!          IIOFF = IIOFF + NSHL(L-1)** 2
!        END IF
!   BUG THE OFFSET SHOULD BE BEFORE THE CYCLE STATEMENT!!
!
! Transform for given L
!
        IF (NSHL(L).GT.0)                                              &
     &  CALL CITRAG(EVEC(1:500),NCF,NVEC,L,NSHL(L),                    &
     &              CIROT(IIOFF),NINL(L),NTESTL,CIOUT, SCR )
!
!   Deallocate storage for this kappa.
!
        CALL DALLOC (JANN, 'JANN', 'MCPIN')
        CALL DALLOC (JBNN, 'JBNN', 'MCPIN')
        CALL DALLOC (INTGRL, 'INTGRL', 'MCPIN')
        CALL DALLOC (CNN, 'CNN', 'MCPIN')
        CALL DALLOC (INTPTR, 'INTPRT', 'MCPIN')

      ENDDO

      CALL dalloc (scr, 'SCR', 'MCPIN')
      CALL dalloc (ciout, 'CIOUT', 'MCPIN')
!
!   Write the rotated CI vectors on file
!
      IF(IBLK.EQ.1) THEN
        J = INDEX(NAME,' ')
        IF (INPCI.EQ.0) THEN
          OPEN (UNIT = 31,FILE=NAME(1:J-1)//'.cbm',FORM='UNFORMATTED',  &
     &    STATUS='UNKNOWN')
        ELSE
          OPEN (UNIT = 31,FILE=NAME(1:J-1)//'.bm',FORM='UNFORMATTED',   &
     &    STATUS='UNKNOWN')
        ENDIF

        WRITE(31) 'G92MIX'
        WRITE(31) NELEC,NCFTOT,NW,NVECTOT,NVECSIZE,NBLOCK
      ENDIF

      WRITE(31) IBLK,NCF,NVEC,IATJPO(1),IASPAR(1)
      WRITE(31) (IVEC(I),I = 1,NVEC)
!     WRITE(31) (IATJPO(I),IASPAR(I),I = 1,NVEC)
      WRITE(31) EAV,(EVAL(I),I = 1,NVEC)
      WRITE(31) ((EVEC(I+(J-1)*NCF),I = 1,NCF),J = 1,NVEC)

      CALL DALLOC (EVAL, 'EVAL', 'MCPIN')
      CALL DALLOC (EVEC, 'EVEC', 'MCPIN')
      CALL DALLOC (IVEC, 'IVEC', 'MCPIN')
      CALL DALLOC (IATJPO, 'IATJPO', 'MCPIN')
      CALL DALLOC (IASPAR, 'IASPAR', 'MCPIN')
 1000 CONTINUE
      CLOSE (NF)
!
!   Close the  mixing  file
!
      CLOSE (30)
      CLOSE (31)
!
      RETURN
      END SUBROUTINE MCPIN

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
Jump to #Navigation Navigation

  City of Davis purchases Boy Scout Cabin, kicks out Boy Scouts
  The Sneeze starts publication.
  Teknolust, a feature film written and directed by UC Davis faculty member Lynn HershmanLeeson is released.
  Friends of the University Airport formed to protect the University Airport.

 May

  Roy Brophy and Alice Gonzales, both former UC Regents, are recipients of the UC Davis Medal on May 20th.

 August

  The Driven to Davis Art in Public Places mural is painted.


text:
!
! Copyright (c) 2005-2018, Christopher C. Hulbert
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
!    list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation
!    and/or other materials provided with the distribution.
!
! 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 HOLDER 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.
!

!   FIXME: Find a less memory hungry way to split the complex data

FUNCTION FMat_VarWriteDouble_0(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(8),INTENT(IN)                       :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteDouble_0

FUNCTION FMat_VarWriteDouble_1(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(8),INTENT(IN),DIMENSION(:)          :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteDouble_1

FUNCTION FMat_VarWriteDouble_2(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(8),INTENT(IN),DIMENSION(:,:)        :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteDouble_2

FUNCTION FMat_VarWriteDouble_3(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(8),INTENT(IN),DIMENSION(:,:,:)      :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteDouble_3

FUNCTION FMat_VarWriteComplexDouble_0(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    DOUBLE COMPLEX,INTENT(IN)                :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteComplexDouble_0

FUNCTION FMat_VarWriteComplexDouble_1(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    DOUBLE COMPLEX,INTENT(IN),DIMENSION(:)   :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress
    DOUBLE PRECISION,DIMENSION(2*SIZE(d))      :: d_split

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    d_split(1:SIZE(d))  = REAL(d)
    d_split(SIZE(d)+1:) = AIMAG(d)
    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d_split,compress)
    END IF
END FUNCTION FMat_VarWriteComplexDouble_1

FUNCTION FMat_VarWriteComplexDouble_2(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    DOUBLE COMPLEX,INTENT(IN),DIMENSION(:,:) :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress
    DOUBLE PRECISION,DIMENSION(2*PRODUCT(SHAPE(d))) :: d_split

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    d_split(1:PRODUCT(SHAPE(d)))  = REAL(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    d_split(PRODUCT(SHAPE(d))+1:) = AIMAG(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d_split,compress)
    END IF
END FUNCTION FMat_VarWriteComplexDouble_2

FUNCTION FMat_VarWriteComplexDouble_3(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                    :: err
    TYPE(mat_t)                                :: mat
    TYPE(matvar_t)                             :: matvar
    DOUBLE COMPLEX,INTENT(IN),DIMENSION(:,:,:) :: d
    INTEGER,INTENT(IN),OPTIONAL                :: compress
    DOUBLE PRECISION,DIMENSION(2*PRODUCT(SHAPE(d))) :: d_split

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    d_split(1:PRODUCT(SHAPE(d)))  = REAL(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    d_split(PRODUCT(SHAPE(d))+1:) = AIMAG(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d_split,compress)
    END IF
END FUNCTION FMat_VarWriteComplexDouble_3

FUNCTION FMat_VarWriteSingle_0(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(4),INTENT(IN)                    :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteSingle_0

FUNCTION FMat_VarWriteSingle_1(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(4),INTENT(IN),DIMENSION(:)          :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteSingle_1

FUNCTION FMat_VarWriteSingle_2(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(4),INTENT(IN),DIMENSION(:,:)        :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteSingle_2

FUNCTION FMat_VarWriteSingle_3(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    REAL(4),INTENT(IN),DIMENSION(:,:,:)      :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteSingle_3

FUNCTION FMat_VarWriteComplexSingle_0(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    COMPLEX(4),INTENT(IN)                    :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteComplexSingle_0

FUNCTION FMat_VarWriteComplexSingle_1(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    COMPLEX(4),INTENT(IN),DIMENSION(:)  :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress
    REAL(4),DIMENSION(2*SIZE(d))             :: d_split

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    d_split(1:SIZE(d))  = REAL(d)
    d_split(SIZE(d)+1:) = AIMAG(d)
    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d_split,compress)
    END IF
END FUNCTION FMat_VarWriteComplexSingle_1

FUNCTION FMat_VarWriteComplexSingle_2(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    COMPLEX(4),INTENT(IN),DIMENSION(:,:) :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress
    REAL(4),DIMENSION(2*PRODUCT(SHAPE(d))) :: d_split

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    d_split(1:PRODUCT(SHAPE(d)))  = REAL(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    d_split(PRODUCT(SHAPE(d))+1:) = AIMAG(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d_split,compress)
    END IF
END FUNCTION FMat_VarWriteComplexSingle_2

FUNCTION FMat_VarWriteComplexSingle_3(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                :: err
    TYPE(mat_t)                            :: mat
    TYPE(matvar_t)                         :: matvar
    COMPLEX(4),INTENT(IN),DIMENSION(:,:,:) :: d
    INTEGER,INTENT(IN),OPTIONAL            :: compress
    REAL(4),DIMENSION(2*PRODUCT(SHAPE(d))) :: d_split

    INTEGER,EXTERNAL                       :: fmat_varwrite_c

    d_split(1:PRODUCT(SHAPE(d)))  = REAL(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    d_split(PRODUCT(SHAPE(d))+1:) = AIMAG(RESHAPE(d,(/2*PRODUCT(SHAPE(d))/)))
    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d_split,compress)
    END IF
END FUNCTION FMat_VarWriteComplexSingle_3

FUNCTION FMat_VarWriteInt32_0(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(4),INTENT(IN)                    :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt32_0

FUNCTION FMat_VarWriteInt32_1(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(4),INTENT(IN),DIMENSION(:)       :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt32_1

FUNCTION FMat_VarWriteInt32_2(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(4),INTENT(IN),DIMENSION(:,:)     :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt32_2

FUNCTION FMat_VarWriteInt32_3(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(4),INTENT(IN),DIMENSION(:,:,:)   :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt32_3

FUNCTION FMat_VarWriteInt16_0(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(2),INTENT(IN)                    :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt16_0

FUNCTION FMat_VarWriteInt16_1(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(2),INTENT(IN),DIMENSION(:)       :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt16_1

FUNCTION FMat_VarWriteInt16_2(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(2),INTENT(IN),DIMENSION(:,:)     :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt16_2

FUNCTION FMat_VarWriteInt16_3(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(2),INTENT(IN),DIMENSION(:,:,:)   :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt16_3

FUNCTION FMat_VarWriteInt8_0(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(1),INTENT(IN)                    :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt8_0

FUNCTION FMat_VarWriteInt8_1(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(1),INTENT(IN),DIMENSION(:)       :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt8_1

FUNCTION FMat_VarWriteInt8_2(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(1),INTENT(IN),DIMENSION(:,:)     :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt8_2

FUNCTION FMat_VarWriteInt8_3(mat,matvar,d,compress) RESULT(err)
IMPLICIT NONE
    INTEGER                                  :: err
    TYPE(mat_t)                              :: mat
    TYPE(matvar_t)                           :: matvar
    INTEGER(1),INTENT(IN),DIMENSION(:,:,:)   :: d
    INTEGER,INTENT(IN),OPTIONAL              :: compress

    INTEGER,EXTERNAL                         :: fmat_varwrite_c

    IF ( .NOT. PRESENT(compress) ) THEN
        err = fmat_varwrite_c(mat,matvar,d,0)
    ELSE
        err = fmat_varwrite_c(mat,matvar,d,compress)
    END IF
END FUNCTION FMat_VarWriteInt8_3

max_stars_repo_name: JohanMabille/matio

max_issues_repo_name: JohanMabille/matio

max_forks_repo_name: JohanMabille/matio


text:
    !> Incorporates (minRestriction + zero average) for fibres network
    !!
	!! @param iShiftDeltaU = nint(CommonPar(1))
	!! @param iShiftLag = nint(CommonPar(2))
	!! @param pen  = commonPar(3)
	!! @param eps  = commonPar(4)

    !! @author Rocha, Felipe Figueredo

!     ------------------------------------------------------------------
Subroutine networkConstraintTheta(AE, BE, MaxLRows, XLL, NDim, iDofT, NodELT, Sol0, &
						Sol1, CommonPar, Param, JParam, DelT, DTm, Time)
!     ------------------------------------------------------------------
	use funcAux , only : numprint, setAEandBE_LM, getSliceAllocate
	use globalVariables, only : NdimE, getMaterial, maxMatPar, getAnisoTensorInv
	use fibresLib
	use fibresMod

	IMPLICIT NONE

	!   ===== SUBROUTINE ARGUMENTS  =======
	integer :: MaxLRows,Ndim,iDofT,NodELT ! all integers
	Real*8 :: DelT, DTm,Time ! all reals
	! Reals Vectors and matrices
	Real*8  :: AE(MaxLRows,MaxLRows), BE(MaxLRows), XLL(*), Sol0(*), Sol1(*), &
		Param(*), JParam(*), CommonPar(*)

	!   =====   END ARGUMENTS  =======

	Integer :: i,j,k,l,m,n,p,q, A, B, ApRow, BpCol , App,Bqp, ApDim, AppDim! counters
	integer, parameter ::  NodG = 2, NodLag = 3 , NdofLag = 6 ! just for 2D
	Real*8 :: MatPar(maxMatPar)
	real*8 , allocatable ::  Xel(:), SolU(:)
	Real*8 :: eps, pen , theta
	Real*8 :: afib(NdimE), Areafib, Lfib, Vf, signal, v(NdimE)
	Real*8 :: MatB(NodG*NdimE,NdofLag), MatC(NdofLag,NodG*NdimE), MatD(NdofLag,NdofLag), VecG(NdofLag)
	integer :: iShiftDeltaU, iShiftU, iShiftLag , iMaterial, constLaw
	integer :: flag1 , flag2, opConstraint

	iShiftU = nint(CommonPar(1))
	iShiftDeltaU = nint(CommonPar(2))
	iShiftLag = nint(CommonPar(3))
	theta = commonPar(4)
	eps = commonPar(5)
	opConstraint = nint(commonPar(6)) ! 0 on boundary and 1 is in all the nodes

	call getSliceAllocate(Xel,XLL,1,NodG,1 ,NdimE, Ndim)
	call get_afib(afib,Lfib,Xel)
	call getSliceAllocate(SolU,Sol1,1,NodG,iShiftU + 1 ,iShiftU + NdimE, iDofT)

	Areafib = Param(Ipos_Areaf)
	afib = Param(Ipos_af:Ipos_af+1)
	Lfib = Param(Ipos_Lf)
	Vf = Param(Ipos_Vf)

	v = -matmul(Bten_invT,afib)/measFibres

!~ 	v = afib

	if(opConstraint == 0) then
		flag1 = nint(Param(Ipos_flag1))
		flag2 = nint(Param(Ipos_flag2))
	else
		flag1 = 1
		flag2 = 1
	end if

	MatC = 0.0d0

	if(flag1 > 0) then
		MatC(1,1) = -v(1)*Areafib
		MatC(2,1) = -v(2)*Areafib
		MatC(3,2) = -v(1)*Areafib
		MatC(4,2) = -v(2)*Areafib
	end if

	if(flag2 > 0) then
		MatC(1,3) = v(1)*Areafib
		MatC(2,3) = v(2)*Areafib
		MatC(3,4) = v(1)*Areafib
		MatC(4,4) = v(2)*Areafib
	end if

	MatC(5,1) = 0.5d0*Vf
	MatC(5,3) = 0.5d0*Vf
	MatC(6,2) = 0.5d0*Vf
	MatC(6,4) = 0.5d0*Vf

	MatC =  theta*MatC
	MatB = transpose(MatC)

	MatD = 0.0d0
	do i = 1,NdofLag
		MatD(i,i) = eps + (1.0d0 - theta)*(1.0d0 - eps)
	end do

	App = (NodLag-1)*iDofT + iShiftLag

	do i = 1, NdofLag
		VecG(i) = (eps + (1.0d0 - theta)*(1.0d0 - eps))*Sol1(App + i)
	end do

	if(iShiftU/=iShiftDeltaU) VecG = VecG - matmul(matC,SolU)

	call setAEandBE_LM(AE,BE,matB,matC,matD,vecG,NodG,NodLag,idofT,NdofLag,iShiftDeltaU,iShiftLag,NdimE)

end subroutine

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Subroutine networkConstraintThetaS(Coupling,CommonPar,iDofT,Ndim,MaxLRows,iAdd)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	use funcAux , only : numprint, setCoupling_lagrangeMultipliers

	use globalVariables, only : NdimE
	IMPLICIT NONE
	Integer :: MaxLRows,Ndim,iDofT, iAdd
	Integer   Coupling(MaxLRows,MaxLRows)
	Real*8 CommonPar(*)
	integer ::  p , q, App, ApRow, BpCol, Bqp, A,B, iShiftDeltaU, iShiftLag , NdofLag
	integer , parameter :: NodG = 2, NodLag = 3

	iShiftDeltaU = nint(CommonPar(2))
	iShiftLag = nint(CommonPar(3))

	NdofLag = NdimE*NdimE + NdimE

	call setCoupling_lagrangeMultipliers(Coupling,NodG,NodLag,idofT,NdofLag,iShiftDeltaU,iShiftLag,NdimE)

end subroutine

max_stars_repo_name: felipefr/gpmaterials

max_issues_repo_name: felipefr/gpmaterials

max_forks_repo_name: felipefr/gpmaterials


text:
	SUBROUTINE G2T_RANGX ( ktype, lunb, nt, itrnd, iret )
C************************************************************************
C* G2T_RANGX								*
C*									*
C* This subroutine adds wind/wave RANGE info from the trending buffer	*
C* to offshore text, OFFTXT.						*
C*									*
C* G2T_RANGX ( KTYPE, LUNB, NT, ITRND, IRET )				*
C*									*
C* Input parameters:							*
C*	KTYPE		INTEGER		Grid parameter type		*
C*					 1 = wave			*
C*					 2 = wind			*
C*	LUNB		INTEGER		LUN for G2T_TXT.TBL		*
C*	NT		INTEGER		Nth time step			*
C*	ITRND		INTEGER		Index for trending buffer	*
C*									*
C* Output Parameters:							*
C*	IRET		INTEGER		Return code			*
C*					 0 = normal return		*
C**									*
C* Log:									*
C* T. Lee/SAIC		03/07						*
C* T. Lee/SAIC		11/07	Added NT for combining periods		*
C* T. Lee/SAIC		06/08	Use 1 FT or less for 0 FT wave		*
C************************************************************************
	INCLUDE		'goftxt.cmn'
	CHARACTER	ooo*72, outstr*72, str(2)*5
	LOGICAL		spread
C-----------------------------------------------------------------------
	iret = 0
C
	DO i = 1, 2
	    str ( i ) = ' '
	END DO
C
	min = mxn_d ( ktype, itrnd, 1 )
	max = mxn_d ( ktype, itrnd, 2 )
	CALL ST_INCH ( min, str ( 1 ), ier )
	CALL ST_INCH ( max, str ( 2 ), ier )
	IF  ( ktype .eq. 1 )  THEN
	    IF ( min .eq. 0 .or. max .eq. 0 )  THEN
		CALL G2T_GTEXT ( lunb, 54, str, ooo, ier )
	      ELSE
		CALL G2T_GTEXT ( lunb, 36, str, ooo, ier )
	    END IF
C
	  ELSE IF ( ktype .eq. 2 )  THEN
	    IF ( max .le. LGHTS )  THEN
		str ( 1 ) = ' '
		CALL ST_INCH ( LGHTS, str ( 2 ), ier )
		CALL G2T_GTEXT ( lunb, 8, str, ooo, ier )
	      ELSE
		IF  ( min .lt. 5 ) THEN
		    str ( 1 ) = ' '
		    CALL G2T_GTEXT ( lunb, 10, str, ooo, ier )
		  ELSE
		    CALL G2T_GTEXT ( lunb, 11, str, ooo, ier )
		END IF
		CALL ST_LSTR ( ooo, looo, ier )
C
		str ( 1 ) = wdir_d ( 1, itrnd )
		str ( 2 ) = wdir_d ( 2, itrnd )
C
C*		The following is the first attempt trying to resolve the
C*		wind spread.  If wind directions spread more than one
C*		cardinal point, use the most prevalent wind direction
C*		only.  Similar code can be found in the G2T_APPEX module.
C
		CALL G2T_WSPREAD ( str, spread, ier )
C
		CALL G2T_GTEXT ( lunb, 9, str, outstr, ier )
		CALL ST_LSTR ( outstr, lout, ier )
		ooo = outstr ( : lout ) // ooo ( : looo )
	    END IF
	END IF
	CALL ST_LSTR ( ooo, looo, ier )
	CALL ST_LSTR ( offtxt ( nt ), lout, ier )
	offtxt ( nt ) = offtxt ( nt ) ( : lout ) // ooo ( : looo )
C
C*      Add location if needed.
C
        CALL ST_LSTR ( id_d ( ktype, itrnd ), lid, ier )
        IF ( idflag .and. ktype .eq. 1 .and. lid .ne. 0 )  THEN
            CALL ST_LSTR ( offtxt ( nt ), loff, ier )
            offtxt ( nt ) = offtxt ( nt ) ( : loff - 1 ) // ' ' //
     +			    id_d ( ktype, itrnd ) ( : lid ) // '.'
        END IF
C*
	RETURN
	END


text:
! RUN: %python %S/test_errors.py %s %flang_fc1
module m1
  implicit none
contains
  subroutine foo(x)
    real :: x
  end subroutine
end module

!Note: PGI, Intel, GNU, and NAG allow this; Sun does not
module m2
  use m1
  implicit none
  !ERROR: 'foo' may not be the name of both a generic interface and a procedure unless it is a specific procedure of the generic
  interface foo
    module procedure s
  end interface
contains
  subroutine s(i)
    integer :: i
  end subroutine
end module

subroutine foo
  !ERROR: Cannot use-associate 'foo'; it is already declared in this scope
  use m1
end

subroutine bar
  !ERROR: Cannot use-associate 'bar'; it is already declared in this scope
  use m1, bar => foo
end

!OK to use-associate a type with the same name as a generic
module m3a
  type :: foo
  end type
end
module m3b
  use m3a
  interface foo
  end interface
end

! Can't have derived type and function with same name
module m4a
  type :: foo
  end type
contains
  !ERROR: 'foo' is already declared in this scoping unit
  function foo(x)
  end
end
! Even if there is also a generic interface of that name
module m4b
  type :: foo
  end type
  !ERROR: 'foo' is already declared in this scoping unit
  interface foo
    procedure :: foo
  end interface foo
contains
  function foo(x)
  end
end
module m4c
  type :: foo
  end type
  interface foo
    !ERROR: 'foo' is already declared in this scoping unit
    real function foo()
    end function foo
  end interface foo
end

! Use associating a name that is a generic and a derived type
module m5a
  interface g
  end interface
  type g
  end type
end module
module m5b
  use m5a
  interface g
    procedure f
  end interface
  type(g) :: x
contains
  function f(i)
  end function
end module
subroutine s5
  use m5b
  type(g) :: y
end

module m6
  real :: f6
  interface g6
  !ERROR: 'f6' is already declared in this scoping unit
    real function f6()
    end function f6
  end interface g6
end module m6

module m7
  integer :: f7
  interface g7
    !ERROR: 'f7' is already declared in this scoping unit
    real function f7()
    end function f7
  end interface g7
end module m7

module m8
  real :: f8
  interface g8
    !ERROR: 'f8' is already declared in this scoping unit
    subroutine f8()
    end subroutine f8
  end interface g8
end module m8

module m9
  type f9
  end type f9
  !ERROR: 'f9' is already declared in this scoping unit
  interface f9
    real function f9()
    end function f9
  end interface f9
contains
  function f9(x)
  end function f9
end module m9

module m10
  type :: t10
  end type t10
  interface f10
    function f10()
    end function f10
  end interface f10
contains
  !ERROR: 'f10' is already declared in this scoping unit
  function f10(x)
  end function f10
end module m10

module m11
  type :: t11
  end type t11
  interface i11
    function f11()
    end function f11
  end interface i11
contains
  !ERROR: 'f11' is already declared in this scoping unit
  function f11(x)
  end function f11
end module m11

module m12
  interface f12
    function f12()
    end function f12
  end interface f12
contains
  !ERROR: 'f12' is already declared in this scoping unit
  function f12(x)
  end function f12
end module m12

module m13
  interface f13
    function f13()
    end function f13
  end interface f13
contains
  !ERROR: 'f13' is already declared in this scoping unit
  function f13()
  end function f13
end module m13

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm