! { dg-do run }
! { dg-options "-std=f2008  -fcheck=bounds" }

! PR fortran/29785
! Check pointer rank remapping at runtime.

! Contributed by Daniel Kraft, d@domob.eu.

PROGRAM main
  IMPLICIT NONE
  INTEGER, TARGET :: arr(12), basem(3, 4)
  INTEGER, POINTER :: vec(:), mat(:, :)
  INTEGER :: i

  arr = (/ (i, i = 1, 12) /)
  basem = RESHAPE (arr, SHAPE (basem))

  ! We need not necessarily change the rank...
  vec(2_1:5) => arr(1_1:12_1:2_1)
  IF (LBOUND (vec, 1) /= 2 .OR. UBOUND (vec, 1) /= 5) STOP 1
  IF (ANY (vec /= (/ 1, 3, 5, 7 /))) STOP 2
  IF (vec(2) /= 1 .OR. vec(5) /= 7) STOP 3

  ! ...but it is of course the more interesting.  Also try remapping a pointer.
  vec => arr(1:12:2)
  mat(1:3, 1:2) => vec
  IF (ANY (LBOUND (mat) /= (/ 1, 1 /) .OR. UBOUND (mat) /= (/ 3, 2 /))) &
    STOP 4
  IF (ANY (mat /= RESHAPE (arr(1:12:2), SHAPE (mat)))) STOP 5
  IF (mat(1, 1) /= 1 .OR. mat(1, 2) /= 7) STOP 6

  ! Remap with target of rank > 1.
  vec(1:12_1) => basem
  IF (LBOUND (vec, 1) /= 1 .OR. UBOUND (vec, 1) /= 12) STOP 7
  IF (ANY (vec /= arr)) STOP 8
  IF (vec(1) /= 1 .OR. vec(5) /= 5 .OR. vec(12) /= 12) STOP 9
END PROGRAM main


text:
C=======================================================================
      PROGRAM COMBTEX
C=======================================================================

C                       Written by Amy Gilkey
C                       Last Revised 04/25/88

C   --COMBTEX combines a TEX file by including (inline) all \include and
C   --\input files.

$HP9000_800 INTRINSICS ON
      CHARACTER*120 INLINE
      CHARACTER*120 FILNAM
      CHARACTER*20 WORD
      LOGICAL INPUT, ABSTRACT, NXTABS

C   --Open the input and output files

      arg_count = iargc()
      if (arg_count .ge. 2) then
         call getarg(2, filnam)
      else
         write (*,*) 'ERROR - No filename specified'
         goto 120
      end if
      NTEX = 20

      OPEN (NTEX, FILE=FILNAM, READONLY, STATUS='OLD', IOSTAT=IERR)
      IF (IERR .NE. 0) THEN
         PRINT *, 'ERROR - Cannot open input file'
         GOTO 120
      END IF

      NCOMB = 19

      OPEN (NCOMB, FILE='OUTFILE', STATUS='NEW')

C   --Process the file

  100 CONTINUE
      ABSTRACT = .FALSE.

      DO WHILE (.TRUE.)
         READ (NTEX, '(Q, A)', END=110) IEND, INLINE
         NCOMB = NCOMB
         INPUT = .FALSE.
         NXTABS = .FALSE.

C      --Scan for \include or \input

         IF (INLINE(1:1) .EQ. '\') THEN
            IE = INDEX (INLINE, '{') - 1
            WORD = INLINE(1:IE)

            IF ((WORD .EQ. '\include') .OR. (WORD .EQ. '\input')
     &         .OR. (WORD .EQ. '\abstract')) THEN
               IF ((WORD .EQ. '\include')
     &            .OR. (WORD .EQ. '\input')) THEN
                  INPUT = .TRUE.
               ELSE IF (WORD .EQ. '\abstract') THEN
                  NXTABS = .TRUE.
               END IF

C            --Change the input file (return to old file at end of file)

               IS = INDEX (INLINE, '{') + 1
               IE = INDEX (INLINE, '}') - 1
               IDIR = INDEX (INLINE, '/') + 1
               I = INDEX (INLINE(IDIR:), '.')
               IF (I .LE. 0) THEN
                  FILNAM = INLINE(IS:IE) // '.tex'
               ELSE
                  FILNAM = INLINE(IS:IE)
               END IF
               IE = INDEX (INLINE, '}')
               PRINT *, INLINE(:IE)

               NTEX = NTEX + 1
               OPEN (NTEX, FILE=FILNAM, READONLY,
     $              STATUS='OLD', IOSTAT=IERR)
               IF (IERR .NE. 0) THEN
                  OPEN (NTEX, FILE=FILNAM//'.tex', READONLY,
     $                 STATUS='OLD', IOSTAT=IERR)
                  IF (IERR .NE. 0) THEN
                     PRINT *, 'ERROR - Cannot open file ', FILNAM
                     NTEX = NTEX - 1
                  end if
               END IF
            END IF
         END IF

C      --Write processed line

         IF (INPUT .OR. ABSTRACT) THEN
            WRITE (NCOMB, 10000) '%%% ', INLINE(:IEND)
         ELSE
            WRITE (NCOMB, 10000) INLINE(:IEND)
         END IF
10000     FORMAT (10A)

         IF (NXTABS) ABSTRACT = .TRUE.
      END DO

  110 CONTINUE

C   --Close TEX file, if input file resume processing with previous file

      IF (NTEX .GT. 20) THEN
         CLOSE (NTEX)
         NTEX = NTEX - 1
         GOTO 100
      END IF

  120 CONTINUE

      CLOSE (NTEX)
      CLOSE (NCOMB)

      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 09.02.03/07                           *
C  *    TEST TITLE : Error indicator = 113                 *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /ERRINF/ ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR,
     1                USRERR, ERRSAV,      FUNSAV,      FILSAV,
     2                EFCNT, EFID
      INTEGER         ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR(10),
     1                USRERR, ERRSAV(200), FUNSAV(200), FILSAV(200),
     2                EFCNT, EFID(100)
      COMMON /ERRCHR/ CURCON,     ERRSRS,    ERRMRK,    ERFLNM,
     1                CONTAB
      CHARACTER       CURCON*200, ERRSRS*40, ERRMRK*20, ERFLNM*80,
     1                CONTAB(40)*150

      INTEGER    IVAL1, SPECWT
      REAL       RA10(10)

      CALL INITGL ('09.02.03/07')

      CALL XPOPPH (ERRFIL, MEMUN)
      CURCON = 'the colour index value is less than zero'
      ERRSRS = '7'
      CALL SETVS ('113', EXPERR, EXPSIZ)

      CALL POPWK (WKID, CONID, WTYPE)
      CALL PQWKC (WKID, ERRIND, CONID, SPECWT)
      CALL CHKINQ ('pqwkc', ERRIND)

      CALL RQCR (WKID, -10, 10, 0, ERRIND, IVAL1, RA10)

      CALL RQPCR (SPECWT, -1, 10, ERRIND, IVAL1, RA10)

      CALL ENDIT

      END


text:
      SUBROUTINE transpose_1(nx,ny,first_x,last_x,input,output)
      integer nx,ny,first_x,last_x
      double complex input(ny,nx),output(nx,ny)
      integer j,k
      !!faster input and slower putput
      DO j=first_x+1,last_x
        DO k=1,ny
          output(j,k)=input(k,j)
        END DO
      END DO
      RETURN
      END

      SUBROUTINE trans_b(nx_th,ny,a,ny_th,nx,b,i0,j0)
      integer nx_th,ny,ny_th,nx
      integer i0,j0
      double complex a(0:ny-1,0:nx_th-1)
      double complex b(0:nx-1,0:ny_th-1)
      integer i,j,ii,jj
      ii=i0
      DO i=0,ny_th-1
        jj=j0
        DO j=0,nx_th-1
        b(jj,i)=a(ii,j)
        jj=jj+1
        END DO
        ii=ii+1
      END DO
      RETURN
      END

      SUBROUTINE transpose_xy(nx,ny,m,j_start,j_end,a,b)
      integer nx,ny,m
      integer j_start, j_end
      double complex a(ny,nx,m), b(nx,ny,m)
      integer i,j,k
      write(*,*) nx, ny, m, j_start,j_end
      DO i=1,m
        DO j=j_start+1,j_end
          DO k=1,nx
            b(k,j,i)=a(j,k,i)
          END DO
        END DO
      END DO
      RETURN
      END

      SUBROUTINE transpose_yx(nx,ny,m,j_start,j_end,a,b)
      integer nx,ny,m
      integer j_start, j_end
      double complex a(nx,ny,m), b(ny,nx,m)
      integer i,j,k
      DO i=1, m
        DO j=j_start+1,j_end
          DO k=1,ny
            b(k,j,i)=a(j,k,i)
          END DO
        END DO
      END DO
      RETURN
      END

max_stars_repo_name: FeiyangJin/hclib

max_issues_repo_name: FeiyangJin/hclib

max_forks_repo_name: FeiyangJin/hclib


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYMAEL
C...Auxiliary to PYSHOW and PYPTFS.
C...Matrix elements for gluon (or photon) emission from
C...a two-body state; to be used by the parton shower routine.
C...Here X_i = 2 E_i/E_cm, R_i = m_i/E_cm and
C...1/sigma_0 d(sigma)/d(x_1)d(x_2) =
C...      = (alpha-strong/2 pi) * CF * PYMAEL,
C...i.e. normalization is such that one recovers the familiar
C...(X1**2+X2**2)/((1-X1)*(1-X2)) for the massless case.
C...Coupling structure:
C...NI =  6- 9 : eikonal soft-gluon expression (spin-independent)
C...   = 11-14 : V -> q qbar (V = vector/axial vector colour singlet)
C...   = 16-19 : q -> q V
C...   = 21-24 : S -> q qbar (S = scalar/pseudoscalar colour singlet)
C...   = 26-29 : q -> q S
C...   = 31-34 : V -> ~q ~qbar  (~q = squark)
C...   = 36-39 : ~q -> ~q V
C...   = 41-44 : S -> ~q ~qbar
C...   = 46-49 : ~q -> ~q S
C...   = 51-54 : chi -> q ~qbar (chi = neutralino/chargino)
C...   = 56-59 : ~q -> q chi
C...   = 61-64 : q -> ~q chi
C...   = 66-69 : ~g -> q ~qbar
C...   = 71-74 : ~q -> q ~g
C...   = 76-79 : q -> ~q ~g
C...   = 81-84 : (9/4)*(eikonal) for gg -> ~g ~g
C...Note that the order of the decay products is important.
C...In each set of four, the variants are ordered as:
C...ICOMBI = 1 : pure non-gamma5, i.e. vector/scalar/...
C...       = 2 : pure gamma5, i.e. axial vector/pseudoscalar/....
C...       = 3 : mixture alpha*(ICOMBI=1) + (1-alpha)*(ICOMBI=2)
C...       = 4 : mixture (ICOMBI=1) +- (ICOMBI=2)

      DOUBLE PRECISION FUNCTION PYMAEL(NI,X1,X2,R1,R2,ALPHA)

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)

C...Check input values. Return zero outside allowed phase space.
      PYMAEL=0D0
      IF(X1.LE.2D0*R1.OR.X1.GE.1D0+R1**2-R2**2) RETURN
      IF(X2.LE.2D0*R2.OR.X2.GE.1D0+R2**2-R1**2) RETURN
      IF(X1+X2.LE.1D0+(R1+R2)**2) RETURN
      IF((2D0-2D0*X1-2D0*X2+X1*X2+2D0*R1**2+2D0*R2**2)**2.GE.
     &(X1**2-4D0*R1**2)*(X2**2-4D0*R2**2)) RETURN
      ALPCOR=MAX(0D0,MIN(1D0,ALPHA))

C...Initial values and flags.
      ICLASS=NI/5
      ICOMBI=NI-5*ICLASS
      ISSET1=0
      ISSET2=0
      ISSET4=0

C... Phase space.
      PS=SQRT((1D0-(R1+R2)**2)*(1D0-(R1-R2)**2))

C...Eikonal expression; also acts as default.
      IF(ICLASS.LE.1.OR.ICLASS.GE.17.OR.ICOMBI.EQ.0) THEN
        RLO=PS
        IF(ICOMBI.EQ.0.OR.ICOMBI.EQ.1) THEN
          ANUM=0D0
        ELSEIF(ICOMBI.EQ.2) THEN
          ANUM=(2D0-X1-X2)**2
        ELSEIF(ICOMBI.EQ.3) THEN
          ANUM=ALPCOR*(2D0-X1-X2)**2
        ELSE
          ANUM=0.5D0*(2D0-X1-X2)**2
        ENDIF
        RFO=PS*2D0*((X1+X2-1D0+ANUM-R1**2-R2**2)/
     &       ((1D0+R1**2-R2**2-X1)*(1D0+R2**2-R1**2-X2))-
     &       R1**2/(1D0+R2**2-R1**2-X2)**2-
     &       R2**2/(1D0+R1**2-R2**2-X1)**2)
        ICOMBI=0

C...V -> q qbar (V = gamma*/Z0/W+-/...).
      ELSEIF(ICLASS.EQ.2) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(2-R1**2-R1**4+6*R1*R2-R2**2+2*R1**2*R2**2-R2**4)/2.D0
        RFO1=-1.D0*(3+6*R1**2+R1**4-6*R1*R2+6*R1**3*R2-2*R2**2
     &       -6*R1**2*R2**2+6*R1*R2**3+R2**4-3*X1+6*R1*R2*X1
     &       +2*R2**2*X1+X1**2-2*R1**2*X1**2+3*R1**2*(2-X1-X2)
     &       +6*R1*R2*(2-X1-X2)-R2**2*(2-X1-X2)-2*X1*(2-X1-X2)
     &       -5*R1**2*X1*(2-X1-X2)+R2**2*X1*(2-X1-X2)+X1**2*(2-X1-X2)
     &       -3*(2-X1-X2)**2-3*R1**2*(2-X1-X2)**2+R2**2*(2-X1-X2)**2
     &       +2*X1*(2-X1-X2)**2+(2-X1-X2)**3-X2)/
     &       (-1+R1**2-R2**2+X2)**2
        RFO1=RFO1-2*(-3+R1**2-6*R1*R2+6*R1**3*R2+3*R2**2-4*R1**2*R2**2
     &       +6*R1*R2**3+2*X1+3*R1**2*X1+R2**2*X1-X1**2-R1**2*X1**2
     &       -R2**2*X1**2+4*(2-X1-X2)+2*R1**2*(2-X1-X2)+3*R1*R2*(2-X1
     &       -X2)-R2**2*(2-X1-X2)-3*X1*(2-X1-X2)-2*R1**2*X1*(2-X1-X2)
     &       +X1**2*(2-X1-X2)-(2-X1-X2)**2-R1**2*(2-X1-X2)**2+R1*R2*(2
     &       -X1-X2)**2+X1*(2-X1-X2)**2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
        RFO1=RFO1-1.D0*(-1+2*R1**2+R1**4+6*R1*R2+6*R1**3*R2-2*R2**2
     &       -6*R1**2*R2**2+6*R1*R2**3+R2**4-X1-2*R1**2*X1-6*R1*R2*X1
     &       +8*R2**2*X1+X1**2-2*R2**2*X1**2-R1**2*(2-X1-X2)+R2**2*(2
     &       -X1-X2)-R1**2*X1*(2-X1-X2)+R2**2*X1*(2-X1-X2)+X1**2*
     &       (2-X1-X2)+X2)/(-1-R1**2+R2**2+X1)**2
        RFO1=RFO1/2.D0
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(2-R1**2-R1**4-6*R1*R2-R2**2+2*R1**2*R2**2-R2**4)/2.D0
        RFO2=-1*(3+6*R1**2+R1**4+6*R1*R2-6*R1**3*R2-2*R2**2
     &       -6*R1**2*R2**2-6*R1*R2**3+R2**4-3*X1-6*R1*R2*X1+2*R2**2*X1
     &       +X1**2-2*R1**2*X1**2+3*R1**2*(2-X1-X2)-6*R1*R2*(2-X1-X2)
     &       -R2**2*(2-X1-X2)-2*X1*(2-X1-X2)-5*R1**2*X1*(2-X1-X2)
     &       +R2**2*X1*(2-X1-X2)+X1**2*(2-X1-X2)-3*(2-X1-X2)**2
     &       -3*R1**2*(2-X1-X2)**2+R2**2*(2-X1-X2)**2+2*X1*(2-X1-X2)**2
     &       +(2-X1-X2)**3-X2)/(-1+R1**2-R2**2+X2)**2
        RFO2=RFO2-2*(-3+R1**2+6*R1*R2-6*R1**3*R2+3*R2**2-4*R1**2*R2**2
     &       -6*R1*R2**3+2*X1+3*R1**2*X1+R2**2*X1-X1**2-R1**2*X1**2
     &       -R2**2*X1**2+4*(2-X1-X2)+2*R1**2*(2-X1-X2)-3*R1*R2*(2-X1
     &       -X2)-R2**2*(2-X1-X2)-3*X1*(2-X1-X2)-2*R1**2*X1*(2-X1-X2)
     &       +X1**2*(2-X1-X2)-(2-X1-X2)**2-R1**2*(2-X1-X2)**2-R1*R2*(2
     &       -X1-X2)**2+X1*(2-X1-X2)**2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
        RFO2=RFO2-1*(-1+2*R1**2+R1**4-6*R1*R2-6*R1**3*R2-2*R2**2
     &       -6*R1**2*R2**2-6*R1*R2**3+R2**4-X1-2*R1**2*X1+6*R1*R2*X1
     &       +8*R2**2*X1+X1**2-2*R2**2*X1**2-R1**2*(2-X1-X2)+R2**2*(2-X1
     &       -X2)-R1**2*X1*(2-X1-X2)+R2**2*X1*(2-X1-X2)+X1**2*(2-X1-X2)
     &       +X2)/(-1-R1**2+R2**2+X1)**2
        RFO2=RFO2/2.D0
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(2D0-R1**2-R1**4-R2**2+2D0*R1**2*R2**2-R2**4)/2D0
        RFO4=(1-R1**4+6*R1**2*R2**2-R2**4+X1+3*R1**2*X1-9*R2**2*X1
     &       -3*X1**2-R1**2*X1**2+3*R2**2*X1**2+X1**3-X2-R1**2*X2
     &       +R2**2*X2-R1**2*X1*X2+R2**2*X1*X2+X1**2*X2)/
     &       (-1-R1**2+R2**2+X1)**2
        RFO4=RFO4
     &       -2*(1+R1**2+R2**2-4*R1**2*R2**2+R1**2*X1+2*R2**2*X1-X1**2
     &       -R2**2*X1**2+2*R1**2*X2+R2**2*X2-3*X1*X2+X1**2*X2-X2**2
     &       -R1**2*X2**2+X1*X2**2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
        RFO4=RFO4+(1-R1**4+6*R1**2*R2**2-R2**4-X1+R1**2*X1-R2**2*X1+X2
     &       -9*R1**2*X2+3*R2**2*X2+R1**2*X1*X2-R2**2*X1*X2-3*X2**2
     &       +3*R1**2*X2**2-R2**2*X2**2+X1*X2**2+X2**3)/
     &       (-1+R1**2-R2**2+X2)**2
        RFO4=RFO4/2.D0
        ISSET4=1
        ENDIF

C...q -> q V.
      ELSEIF(ICLASS.EQ.3) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0-2D0*R1**2+R1**4+R2**2-6D0*R1*R2**2
     &        +R1**2*R2**2-2D0*R2**4)
        RFO1=2*(-1+R1-2*R1**2+2*R1**3-R1**4+R1**5-R2**2+R1*R2**2
     &       -5*R1**2*R2**2+R1**3*R2**2-2*R1*R2**4+2*X1-2*R1*X1
     &       +2*R1**2*X1-2*R1**3*X1+2*R2**2*X1+5*R1*R2**2*X1
     &       +R1**2*R2**2*X1+2*R2**4*X1-X1**2+R1*X1**2-R2**2*X1**2+3*X2
     &       +4*R1**2*X2+R1**4*X2+2*R2**2*X2+2*R1**2*R2**2*X2-4*X1*X2
     &       -2*R1**2*X1*X2-R2**2*X1*X2+X1**2*X2-2*X2**2
     &       -2*R1**2*X2**2+X1*X2**2)/(1-R1**2+R2**2-X2)/(-2+X1+X2)
        RFO1=RFO1+(2*R2**2+6*R1*R2**2-6*R1**2*R2**2+6*R1**3*R2**2
     &       +2*R2**4+6*R1*R2**4-R2**2*X1+R1**2*R2**2*X1-R2**4*X1+X2
     &       -R1**4*X2-3*R2**2*X2-6*R1*R2**2*X2+9*R1**2*R2**2*X2
     &       -2*R2**4*X2-X1*X2+R1**2*X1*X2-X2**2-3*R1**2*X2**2
     &       +2*R2**2*X2**2+X1*X2**2)/(-1+R1**2-R2**2+X2)**2
        RFO1=RFO1+(-4-8*R1**2-4*R1**4+4*R2**2-4*R1**2*R2**2+8*R2**4
     &       +9*X1+10*R1**2*X1+R1**4*X1-3*R2**2*X1+6*R1*R2**2*X1
     &       +R1**2*R2**2*X1-2*R2**4*X1-6*X1**2-2*R1**2*X1**2+X1**3
     &       +7*X2+8*R1**2*X2+R1**4*X2-7*R2**2*X2+6*R1*R2**2*X2
     &       +R1**2*R2**2*X2-2*R2**4*X2-9*X1*X2-3*R1**2*X1*X2
     &       +2*R2**2*X1*X2+2*X1**2*X2-3*X2**2-R1**2*X2**2
     &       +2*R2**2*X2**2+X1*X2**2)/(-2+X1+X2)**2
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0-2D0*R1**2+R1**4+R2**2+6D0*R1*R2**2
     &        +R1**2*R2**2-2D0*R2**4)
        RFO2=2*(1+R1+2*R1**2+2*R1**3+R1**4+R1**5+R2**2+R1*R2**2
     &       +5*R1**2*R2**2+R1**3*R2**2-2*R1*R2**4-2*X1-2*R1*X1
     &       -2*R1**2*X1-2*R1**3*X1-2*R2**2*X1+5*R1*R2**2*X1
     &       -R1**2*R2**2*X1-2*R2**4*X1+X1**2+R1*X1**2+R2**2*X1**2-3*X2
     &       -4*R1**2*X2-R1**4*X2-2*R2**2*X2-2*R1**2*R2**2*X2+4*X1*X2
     &       +2*R1**2*X1*X2+R2**2*X1*X2-X1**2*X2+2*X2**2+2*R1**2*X2**2
     &       -X1*X2**2)/(-1+R1**2-R2**2+X2)/(-2+X1+X2)
        RFO2=RFO2+(2*R2**2-6*R1*R2**2-6*R1**2*R2**2-6*R1**3*R2**2
     &       +2*R2**4-6*R1*R2**4-R2**2*X1+R1**2*R2**2*X1-R2**4*X1+X2
     &       -R1**4*X2-3*R2**2*X2+6*R1*R2**2*X2+9*R1**2*R2**2*X2
     &       -2*R2**4*X2-X1*X2+R1**2*X1*X2-X2**2-3*R1**2*X2**2
     &       +2*R2**2*X2**2+X1*X2**2)/(-1+R1**2-R2**2+X2)**2
        RFO2=RFO2+(-4-8*R1**2-4*R1**4+4*R2**2-4*R1**2*R2**2+8*R2**4+9*X1
     &       +10*R1**2*X1+R1**4*X1-3*R2**2*X1-6*R1*R2**2*X1
     &       +R1**2*R2**2*X1-2*R2**4*X1-6*X1**2-2*R1**2*X1**2+X1**3
     &       +7*X2+8*R1**2*X2+R1**4*X2-7*R2**2*X2-6*R1*R2**2*X2
     &       +R1**2*R2**2*X2-2*R2**4*X2-9*X1*X2-3*R1**2*X1*X2
     &       +2*R2**2*X1*X2+2*X1**2*X2-3*X2**2-R1**2*X2**2+2*R2**2*X2**2
     &       +X1*X2**2)/(-2+X1+X2)**2
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1.D0-2.D0*R1**2+R1**4+R2**2+R1**2*R2**2-2.D0*R2**4)
        RFO4=2*(1+2*R1**2+R1**4+R2**2+5*R1**2*R2**2-2*X1-2*R1**2*X1
     &       -2*R2**2*X1-R1**2*R2**2*X1-2*R2**4*X1+X1**2+R2**2*X1**2
     &       -3*X2-4*R1**2*X2-R1**4*X2-2*R2**2*X2-2*R1**2*R2**2*X2
     &       +4*X1*X2+2*R1**2*X1*X2+R2**2*X1*X2-X1**2*X2+2*X2**2
     &       +2*R1**2*X2**2-X1*X2**2)/(-1+R1**2-R2**2+X2)/(-2+X1+X2)
        RFO4=RFO4+(2*R2**2-6*R1**2*R2**2+2*R2**4-R2**2*X1+R1**2*R2**2*X1
     &       -R2**4*X1+X2-R1**4*X2-3*R2**2*X2+9*R1**2*R2**2*X2
     &       -2*R2**4*X2-X1*X2+R1**2*X1*X2-X2**2-3*R1**2*X2**2
     &       +2*R2**2*X2**2+X1*X2**2)/(-1+R1**2-R2**2+X2)**2
        RFO4=RFO4+(-4-8*R1**2-4*R1**4+4*R2**2-4*R1**2*R2**2+8*R2**4+9*X1
     &       +10*R1**2*X1+R1**4*X1-3*R2**2*X1+R1**2*R2**2*X1-2*R2**4*X1
     &       -6*X1**2-2*R1**2*X1**2+X1**3+7*X2+8*R1**2*X2+R1**4*X2
     &       -7*R2**2*X2+R1**2*R2**2*X2-2*R2**4*X2-9*X1*X2-3*R1**2*X1*X2
     &       +2*R2**2*X1*X2+2*X1**2*X2-3*X2**2-R1**2*X2**2+2*R2**2*X2**2
     &       +X1*X2**2)/(2-X1-X2)**2
        ISSET4=1
        ENDIF

C...S -> q qbar    (S = h0/H0/A0/H+-/...).
      ELSEIF(ICLASS.EQ.4) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0-R1**2-R2**2-2D0*R1*R2)
        RFO1=-(-1+R1**4-2*R1*R2-2*R1**3*R2-6*R1**2*R2**2-2*R1*R2**3
     &       +R2**4+X1-R1**2*X1+2*R1*R2*X1+3*R2**2*X1+X2+R1**2*X2
     &       -R2**2*X2-X1*X2)/(-1-R1**2+R2**2+X1)**2
     &       -2*(R1**2+R1**4-2*R1**3*R2+R2**2-6*R1**2*R2**2-2*R1*R2**3
     &       +R2**4-R1**2*X1+R1*R2*X1+2*R2**2*X1+2*R1**2*X2+R1*R2*X2
     &       -R2**2*X2-X1*X2)/(-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
     &       -(-1+R1**4-2*R1*R2-2*R1**3*R2-6*R1**2*R2**2-2*R1*R2**3
     &       +R2**4+X1-R1**2*X1+R2**2*X1+X2+3*R1**2*X2+2*R1*R2*X2
     &       -R2**2*X2-X1*X2)/(-1+R1**2-R2**2+X2)**2
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0-R1**2-R2**2+2D0*R1*R2)
        RFO2=-(-1+R1**4+2*R1*R2+2*R1**3*R2-6*R1**2*R2**2+2*R1*R2**3
     &       +R2**4+X1-R1**2*X1-2*R1*R2*X1+3*R2**2*X1+X2+R1**2*X2
     &       -R2**2*X2-X1*X2)/(-1-R1**2+R2**2+X1)**2
     &       -(-1+R1**4+2*R1*R2+2*R1**3*R2-6*R1**2*R2**2+2*R1*R2**3
     &       +R2**4+X1-R1**2*X1+R2**2*X1+X2+3*R1**2*X2-2*R1*R2*X2
     &       -R2**2*X2-X1*X2)/(-1+R1**2-R2**2+X2)**2
     &       +2*(-R1**2-R1**4-2*R1**3*R2-R2**2+6*R1**2*R2**2
     &       -2*R1*R2**3-R2**4+R1**2*X1+R1*R2*X1-2*R2**2*X1
     &       -2*R1**2*X2+R1*R2*X2+R2**2*X2+X1*X2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1D0-R1**2-R2**2)
        RFO4=-(-1+R1**4-6*R1**2*R2**2+R2**4+X1-R1**2*X1+3*R2**2*X1+X2
     &       +R1**2*X2-R2**2*X2-X1*X2)/(-1-R1**2+R2**2+X1)**2
     &       -2*(R1**2+R1**4+R2**2-6*R1**2*R2**2+R2**4-R1**2*X1
     &       +2*R2**2*X1+2*R1**2*X2-R2**2*X2-X1*X2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
     &       -(-1+R1**4-6*R1**2*R2**2+R2**4+X1-R1**2*X1+R2**2*X1
     &       +X2+3*R1**2*X2-R2**2*X2-X1*X2)/(-1+R1**2-R2**2+X2)**2
        ISSET4=1
        ENDIF

C...q -> q S.
      ELSEIF(ICLASS.EQ.5) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0+R1**2-R2**2+2D0*R1)
        RFO1=(4-4*R1**2+4*R2**2-3*X1-2*R1*X1+R1**2*X1-R2**2*X1-5*X2
     &       -2*R1*X2+R1**2*X2-R2**2*X2+X1*X2+X2**2)/(-2+X1+X2)**2
     &       +2*(3-R1-5*R1**2-R1**3+3*R2**2+R1*R2**2-2*X1-R1*X1
     &       +R1**2*X1-4*X2+2*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (1-R1**2+R2**2-X2)/(-2+X1+X2)
     &       +(2-2*R1-6*R1**2-2*R1**3+2*R2**2-2*R1*R2**2-X1+R1**2*X1
     &       -R2**2*X1-3*X2+2*R1*X2+3*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (-1+R1**2-R2**2+X2)**2
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0+R1**2-R2**2-2D0*R1)
        RFO2=(4-4*R1**2+4*R2**2-3*X1+2*R1*X1+R1**2*X1-R2**2*X1-5*X2
     &       +2*R1*X2+R1**2*X2-R2**2*X2+X1*X2+X2**2)/(-2+X1+X2)**2
     &       +2*(3+R1-5*R1**2+R1**3+3*R2**2-R1*R2**2-2*X1+R1*X1
     &       +R1**2*X1-4*X2+2*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (1-R1**2+R2**2-X2)/(-2+X1+X2)
     &       +(2+2*R1-6*R1**2+2*R1**3+2*R2**2+2*R1*R2**2-X1+R1**2*X1
     &       -R2**2*X1-3*X2-2*R1*X2+3*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (-1+R1**2-R2**2+X2)**2
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1D0+R1**2-R2**2)
        RFO4=(4-4*R1**2+4*R2**2-3*X1+R1**2*X1-R2**2*X1-5*X2+R1**2*X2
     &       -R2**2*X2+X1*X2+X2**2)/(-2+X1+X2)**2
     &       +2*(3-5*R1**2+3*R2**2-2*X1+R1**2*X1-4*X2+2*R1**2*X2
     &       -R2**2*X2+X1*X2+X2**2)/(1-R1**2+R2**2-X2)/(-2+X1+X2)
     &       +(2-6*R1**2+2*R2**2-X1+R1**2*X1-R2**2*X1-3*X2+3*R1**2*X2
     &       -R2**2*X2+X1*X2+X2**2)/(-1+R1**2-R2**2+X2)**2
        ISSET4=1
        ENDIF

C...V -> ~q ~qbar  (~q = squark).
      ELSEIF(ICLASS.EQ.6) THEN
        RLO1=PS*(1D0-2D0*R1**2+R1**4-2D0*R2**2-2D0*R1**2*R2**2+R2**4)
        RFO1=2D0*3D0+(1+R1**2+R2**2-X1)*(4*R1**2-X1**2)/
     &       (-1-R1**2+R2**2+X1)**2
     &       -2D0*(-1-3*R1**2-R2**2+X1+X1**2/2+X2-X1*X2/2)/
     &       (-1-R1**2+R2**2+X1)
     &       +(1+R1**2+R2**2-X2)*(4*R2**2-X2**2)
     &       /(-1+R1**2-R2**2+X2)**2
     &       -2D0*(-1-R1**2-3*R2**2+X1+X2-X1*X2/2+X2**2/2)/
     &       (-1+R1**2-R2**2+X2)
     &       -(-4*R1**2-4*R1**4-4*R2**2-8*R1**2*R2**2-4*R2**4+2*X1
     &       +6*R1**2*X1+6*R2**2*X1-2*X1**2+2*X2+6*R1**2*X2+6*R2**2*X2
     &       -4*X1*X2-2*R1**2*X1*X2-2*R2**2*X1*X2+X1**2*X2-2*X2**2
     &       +X1*X2**2)/(-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
        ISSET1=1

C...~q -> ~q V.
      ELSEIF(ICLASS.EQ.7) THEN
        RLO1=PS*(1D0-2D0*R1**2+R1**4-2D0*R2**2-2D0*R1**2*R2**2+R2**4)
        RFO1=16*R2**2+8*(4*R2**2+2*R2**2*X1+X2+R1**2*X2+R2**2*X2-X1*X2
     &       -2*X2**2)/(3*(-1+R1**2-R2**2+X2))+8*(1+R1**2+R2**2-X2)*
     &       (4*R2**2-X2**2)/(3*(-1+R1**2-R2**2+X2)**2)+8*(X1+X2)*
     &       (-1-2*R1**2-R1**4-2*R2**2+2*R1**2*R2**2-R2**4+2*X1
     &       +2*R1**2*X1+2*R2**2*X1-X1**2+2*X2+2*R1**2*X2+2*R2**2*X2
     &       -2*X1*X2-X2**2)/(3*(-2+X1+X2)**2)+8*(-1-R1**2+R2**2-X1)*
     &       (2*R2**2*X1+X2+R1**2*X2+R2**2*X2-X1*X2-X2**2)/
     &       (3*(-1+R1**2-R2**2+X2)*(-2+X1+X2))+8*(1+2*R1**2+R1**4
     &       +2*R2**2-2*R1**2*R2**2+R2**4-2*X1-2*R1**2*X1-4*R2**2*X1
     &       +X1**2-3*X2-3*R1**2*X2-3*R2**2*X2+3*X1*X2+2*X2**2)/
     &       (3*(-2+X1+X2))
        RFO1=3D0*RFO1/8D0
        ISSET1=1

C...S -> ~q ~qbar.
      ELSEIF(ICLASS.EQ.8) THEN
        RLO1=PS
        RFO1=(-1-2*R1**2-R1**4-2*R2**2+2*R1**2*R2**2-R2**4+2*X1
     &       +2*R1**2*X1+2*R2**2*X1-X1**2-R2**2*X1**2+2*X2+2*R1**2*X2
     &       +2*R2**2*X2-3*X1*X2-R1**2*X1*X2-R2**2*X1*X2+X1**2*X2-X2**2
     &       -R1**2*X2**2+X1*X2**2)/
     &       (1+R1**2-R2**2-X1)**2/(-1+R1**2-R2**2+X2)**2
        RFO1=2D0*RFO1
        ISSET1=1

C...~q -> ~q S.
      ELSEIF(ICLASS.EQ.9) THEN
        RLO1=PS
        RFO1=(-1-R1**2-R2**2+X2)/(-1+R1**2-R2**2+X2)**2
     &       +(1+R1**2-R2**2+X1)/(-1+R1**2-R2**2+X2)/(-2+X1+X2)
     &       -(X1+X2)/(-2+X1+X2)**2
        ISSET1=1

C...chi -> q ~qbar   (chi = neutralino/chargino).
      ELSEIF(ICLASS.EQ.10) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0+R1**2-R2**2+2D0*R1)
        RFO1=(2*R1+X1)*(-1-R1**2-R2**2+X1)/(-1-R1**2+R2**2+X1)**2
     &       +2*(-1-R1**2-2*R1**3-R2**2-2*R1*R2**2+3*X1/2+R1*X1
     &       -R1**2*X1/2-R2**2*X1/2+X2+R1*X2+R1**2*X2-X1*X2/2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
     &       +(2-2*R1-6*R1**2-2*R1**3+2*R2**2-2*R1*R2**2-X1+R1**2*X1
     &       -R2**2*X1-3*X2+2*R1*X2+3*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (-1+R1**2-R2**2+X2)**2
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0-2D0*R1+R1**2-R2**2)
        RFO2=(2*R1-X1)*(1+R1**2+R2**2-X1)/(-1-R1**2+R2**2+X1)**2
     &       +2*(-1-R1**2+2*R1**3-R2**2+2*R1*R2**2+3*X1/2-R1*X1
     &       -R1**2*X1/2-R2**2*X1/2+X2-R1*X2+R1**2*X2-X1*X2/2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
     &       +(2+2*R1-6*R1**2+2*R1**3+2*R2**2+2*R1*R2**2-X1+R1**2*X1
     &       -R2**2*X1-3*X2-2*R1*X2+3*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (-1+R1**2-R2**2+X2)**2
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1+R1**2-R2**2)
        RFO4=X1*(-1-R1**2-R2**2+X1)/(-1-R1**2+R2**2+X1)**2
     &       +2D0*(-1-R1**2-R2**2+3*X1/2-R1**2*X1/2-R2**2*X1/2
     &       +X2+R1**2*X2-X1*X2/2)/
     &       (-1-R1**2+R2**2+X1)/(-1+R1**2-R2**2+X2)
     &       +(2-6*R1**2+2*R2**2-X1+R1**2*X1-R2**2*X1-3*X2+3*R1**2*X2
     &       -R2**2*X2+X1*X2+X2**2)/(-1+R1**2-R2**2+X2)**2
        ISSET4=1
        ENDIF

C...~q -> q chi.
      ELSEIF(ICLASS.EQ.11) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0-(R1+R2)**2)
        RFO1=(1+R1**2+2*R1*R2+R2**2-X1-X2)*(X1+X2)/(-2+X1+X2)**2
     &       -(-1+R1**4-2*R1*R2-2*R1**3*R2-6*R1**2*R2**2-2*R1*R2**3
     &       +R2**4+X1-R1**2*X1+R2**2*X1+X2+3*R1**2*X2+2*R1*R2*X2
     &       -R2**2*X2-X1*X2)/(-1+R1**2-R2**2+X2)**2
     &       +(-1-2*R1**2-R1**4-2*R1*R2-2*R1**3*R2+2*R1*R2**3+R2**4
     &       +X1+R1**2*X1-2*R1*R2*X1-3*R2**2*X1+2*R1**2*X2-2*R2**2*X2
     &       +X1*X2+X2**2)/(-1+R1**2-R2**2+X2)/(-2+X1+X2)
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0-(R1-R2)**2)
        RFO2=(1+R1**2-2*R1*R2+R2**2-X1-X2)*(X1+X2)/
     &       (-2+X1+X2)**2
     &       -(-1+R1**4+2*R1*R2+2*R1**3*R2-6*R1**2*R2**2+2*R1*R2**3
     &       +R2**4+X1-R1**2*X1+R2**2*X1+X2+3*R1**2*X2-2*R1*R2*X2
     &       -R2**2*X2-X1*X2)/(-1+R1**2-R2**2+X2)**2
     &       +(-1-2*R1**2-R1**4+2*R1*R2+2*R1**3*R2-2*R1*R2**3+R2**4
     &       +X1+R1**2*X1+2*R1*R2*X1-3*R2**2*X1+2*R1**2*X2-2*R2**2*X2
     &       +X1*X2+X2**2)/(-1+R1**2-R2**2+X2)/(-2+X1+X2)
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1D0-R1**2-R2**2)
        RFO4=(1+R1**2+R2**2-X1-X2)*(X1+X2)/(-2+X1+X2)**2
     &       -(-1+R1**4-6*R1**2*R2**2+R2**4+X1-R1**2*X1+R2**2*X1+X2
     &       +3*R1**2*X2-R2**2*X2-X1*X2)/
     &       (-1+R1**2-R2**2+X2)**2
     &       -(-1-2*R1**2-R1**4+R2**4+X1+R1**2*X1-3*R2**2*X1
     &       +2*R1**2*X2-2*R2**2*X2+X1*X2+X2**2)/
     &       (2-X1-X2)/(-1+R1**2-R2**2+X2)
        ISSET4=1
        ENDIF

C...q -> ~q chi.
      ELSEIF(ICLASS.EQ.12) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0-R1**2+R2**2+2D0*R2)
        RFO1=(2*R2+X2)*(-1-R1**2-R2**2+X2)/(-1+R1**2-R2**2+X2)**2
     &       +(4+4*R1**2-4*R2**2-5*X1-R1**2*X1-2*R2*X1+R2**2*X1+X1**2
     &       -3*X2-R1**2*X2-2*R2*X2+R2**2*X2+X1*X2)/
     &       (-2+X1+X2)**2-2*(-1-R1**2+R2+R1**2*R2-R2**2-R2**3+X1
     &       +R2*X1+R2**2*X1+2*X2+R1**2*X2-X1*X2/2-X2**2/2)/
     &       (2-X1-X2)/(-1+R1**2-R2**2+X2)
        ISSET1=1
        END IF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0-R1**2+R2**2-2D0*R2)
        RFO2=(2*R2-X2)*(1+R1**2+R2**2-X2)/(-1+R1**2-R2**2+X2)**2
     &       +(4+4*R1**2-4*R2**2-5*X1-R1**2*X1+2*R2*X1+R2**2*X1+X1**2
     &       -3*X2-R1**2*X2+2*R2*X2+R2**2*X2+X1*X2)/
     &       (-2+X1+X2)**2-2*(-1-R1**2-R2-R1**2*R2-R2**2+R2**3+X1
     &       -R2*X1+R2**2*X1+2*X2+R1**2*X2-X1*X2/2-X2**2/2)/
     &       (2-X1-X2)/(-1+R1**2-R2**2+X2)
        ISSET2=1
        END IF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1D0-R1**2+R2**2)
        RFO4=X2*(-1-R1**2-R2**2+X2)/(-1+R1**2-R2**2+X2)**2
     &       +(4+4*R1**2-4*R2**2-5*X1-R1**2*X1+R2**2*X1+X1**2
     &       -3*X2-R1**2*X2+R2**2*X2+X1*X2)/
     &       (-2+X1+X2)**2-2*(-1-R1**2-R2**2+X1+R2**2*X1+2*X2
     &       +R1**2*X2-X1*X2/2-X2**2/2)/
     &       (2-X1-X2)/(-1+R1**2-R2**2+X2)
        ISSET4=1
        END IF

C...~g -> q ~qbar.
      ELSEIF(ICLASS.EQ.13) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0+R1**2-R2**2+2D0*R1)
        RFO1=4*(2*R1+X1)*(-1-R1**2-R2**2+X1)/(3*(-1-R1**2+R2**2+X1)**2)
     &       -(-1-R1**2-2*R1**3-R2**2-2*R1*R2**2+3*X1/2+R1*X1-R1**2*X1/2
     &       -R2**2*X1/2+X2+R1*X2+R1**2*X2-X1*X2/2)/(3*(-1-R1**2+R2**2
     &       +X1)*(-1+R1**2-R2**2+X2))-3*(-1+R1-R1**2-R1**3-R2**2
     &       +R1*R2**2+2*X1+R2**2*X1-X1**2/2+X2+R1*X2+R1**2*X2-X1*X2/2)/
     &       ((-1-R1**2+R2**2+X1)*(2-X1-X2))+3*(4-4*R1**2+4*R2**2-3*X1
     &       -2*R1*X1+R1**2*X1-R2**2*X1-5*X2-2*R1*X2+R1**2*X2-R2**2*X2
     &       +X1*X2+X2**2)/(-2+X1+X2)**2+3*(3-R1-5*R1**2-R1**3+3*R2**2
     &       +R1*R2**2-2*X1-R1*X1+R1**2*X1-4*X2+2*R1**2*X2-R2**2*X2
     &       +X1*X2+X2**2)/((1-R1**2+R2**2-X2)*(-2+X1+X2))+4*(2-2*R1
     &       -6*R1**2-2*R1**3+2*R2**2-2*R1*R2**2-X1+R1**2*X1-R2**2*X1
     &       -3*X2+2*R1*X2+3*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (3*(-1+R1**2-R2**2+X2)**2)
        RFO1=3D0*RFO1/4D0
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0+R1**2-R2**2-2D0*R1)
        RFO2=4*(2*R1-X1)*(1+R1**2+R2**2-X1)/(3*(-1-R1**2+R2**2+X1)**2)
     &       -3*(-1-R1-R1**2+R1**3-R2**2-R1*R2**2+2*X1+R2**2*X1-X1**2/2
     &       +X2-R1*X2+R1**2*X2-X1*X2/2)/((-1-R1**2+R2**2+X1)*(2-X1-X2))
     &       +(2+2*R1**2-4*R1**3+2*R2**2-4*R1*R2**2-3*X1+2*R1*X1
     &       +R1**2*X1+R2**2*X1-2*X2+2*R1*X2-2*R1**2*X2+X1*X2)/
     &       (6*(-1-R1**2+R2**2+X1)*(-1+R1**2-R2**2+X2))+3*(4-4*R1**2
     &       +4*R2**2-3*X1+2*R1*X1+R1**2*X1-R2**2*X1-5*X2+2*R1*X2
     &       +R1**2*X2-R2**2*X2+X1*X2+X2**2)/(-2+X1+X2)**2+3*(3+R1
     &       -5*R1**2+R1**3+3*R2**2-R1*R2**2-2*X1+R1*X1+R1**2*X1-4*X2
     &       +2*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       ((1-R1**2+R2**2-X2)*(-2+X1+X2))+4*(2+2*R1-6*R1**2+2*R1**3
     &       +2*R2**2+2*R1*R2**2-X1+R1**2*X1-R2**2*X1-3*X2-2*R1*X2
     &       +3*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (3*(-1+R1**2-R2**2+X2)**2)
        RFO2=3D0*RFO2/4D0
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1D0+R1**2-R2**2)
        RFO4=8*X1*(-1-R1**2-R2**2+X1)/(3*(-1-R1**2+R2**2+X1)**2)-6*(-1
     &       -R1**2-R2**2+2*X1+R2**2*X1-X1**2/2+X2+R1**2*X2-X1*X2/2)/
     &       ((-1-R1**2+R2**2+X1)*(2-X1-X2))+(2+2*R1**2+2*R2**2-3*X1
     &       +R1**2*X1+R2**2*X1-2*X2-2*R1**2*X2+X1*X2)/(3*(-1-R1**2
     &       +R2**2+X1)*(-1+R1**2-R2**2+X2))+6*(4-4*R1**2+4*R2**2-3*X1
     &       +R1**2*X1-R2**2*X1-5*X2+R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (-2+X1+X2)**2+6*(3-5*R1**2+3*R2**2-2*X1+R1**2*X1-4*X2
     &       +2*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       ((1-R1**2+R2**2-X2)*(-2+X1+X2))+8*(2-6*R1**2+2*R2**2-X1
     &       +R1**2*X1-R2**2*X1-3*X2+3*R1**2*X2-R2**2*X2+X1*X2+X2**2)/
     &       (3*(-1+R1**2-R2**2+X2)**2)
        RFO4=3D0*RFO4/8D0
        ISSET4=1
        ENDIF

C...~q -> q ~g.
      ELSEIF(ICLASS.EQ.14) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1-R1**2-R2**2-2D0*R1*R2)
        RFO1=64*(1+R1**2+2*R1*R2+R2**2-X1-X2)*(X1+X2)/(9*(-2+X1+X2)**2)
     &       -16*(-1+R1**4-2*R1*R2-2*R1**3*R2-6*R1**2*R2**2-2*R1*R2**3
     &       +R2**4+X1-R1**2*X1+2*R1*R2*X1+3*R2**2*X1+X2+R1**2*X2
     &       -R2**2*X2-X1*X2)/(-1-R1**2+R2**2+X1)**2-16*(R1**2+R1**4
     &       -2*R1**3*R2+R2**2-6*R1**2*R2**2-2*R1*R2**3+R2**4
     &       -R1**2*X1+R1*R2*X1+2*R2**2*X1+2*R1**2*X2+R1*R2*X2-R2**2*X2
     &       -X1*X2)/((-1-R1**2+R2**2+X1)*(-1+R1**2-R2**2+X2))
     &       -64*(-1+R1**4-2*R1*R2-2*R1**3*R2-6*R1**2*R2**2-2*R1*R2**3
     &       +R2**4+X1-R1**2*X1+R2**2*X1+X2+3*R1**2*X2+2*R1*R2*X2
     &       -R2**2*X2-X1*X2)/(9*(-1+R1**2-R2**2+X2)**2)
     &       +8*(-1+R1**4-2*R1*R2+2*R1**3*R2-2*R2**2-2*R1*R2**3-R2**4
     &       -2*R1**2*X1+2*R2**2*X1+X1**2+X2-3*R1**2*X2-2*R1*R2*X2
     &       +R2**2*X2+X1*X2)/((-1-R1**2+R2**2+X1)*(-2+X1+X2))
        RFO1=RFO1
     &       +8*(-1-2*R1**2-R1**4-2*R1*R2-2*R1**3*R2+2*R1*R2**3+R2**4
     &       +X1+R1**2*X1-2*R1*R2*X1-3*R2**2*X1+2*R1**2*X2-2*R2**2*X2
     &       +X1*X2+X2**2)/(9*(2-X1-X2)*(-1+R1**2-R2**2+X2))
        RFO1=9D0*RFO1/64D0
        ISSET1=1
        ENDIF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1-R1**2-R2**2+2D0*R1*R2)
        RFO2=64*(1+R1**2-2*R1*R2+R2**2-X1-X2)*(X1+X2)/(9*(-2+X1+X2)**2)
     &       -16*(-1+R1**4+2*R1*R2+2*R1**3*R2-6*R1**2*R2**2+2*R1*R2**3
     &       +R2**4+X1-R1**2*X1-2*R1*R2*X1+3*R2**2*X1+X2+R1**2*X2
     &       -R2**2*X2-X1*X2)/(-1-R1**2+R2**2+X1)**2-64*(-1+R1**4
     &       +2*R1*R2+2*R1**3*R2-6*R1**2*R2**2+2*R1*R2**3+R2**4+X1
     &       -R1**2*X1+R2**2*X1+X2+3*R1**2*X2-2*R1*R2*X2-R2**2*X2
     &       -X1*X2)/(9*(-1+R1**2-R2**2+X2)**2)+16*(-R1**2-R1**4
     &       -2*R1**3*R2-R2**2+6*R1**2*R2**2-2*R1*R2**3-R2**4+R1**2*X1
     &       +R1*R2*X1-2*R2**2*X1-2*R1**2*X2+R1*R2*X2+R2**2*X2+X1*X2)/
     &       ((-1-R1**2+R2**2+X1)*(-1+R1**2-R2**2+X2))
        RFO2=RFO2
     &       +8*(-1+R1**4+2*R1*R2-2*R1**3*R2-2*R2**2+2*R1*R2**3-R2**4
     &       -2*R1**2*X1+2*R2**2*X1+X1**2+X2-3*R1**2*X2+2*R1*R2*X2
     &       +R2**2*X2+X1*X2)/((-1-R1**2+R2**2+X1)*(-2+X1+X2))
     &       +8*(-1-2*R1**2-R1**4+2*R1*R2+2*R1**3*R2-2*R1*R2**3
     &       +R2**4+X1+R1**2*X1+2*R1*R2*X1-3*R2**2*X1+2*R1**2*X2
     &       -2*R2**2*X2+X1*X2+X2**2)/(9*(2-X1-X2)*(-1+R1**2-R2**2+X2))
        RFO2=9D0*RFO2/64D0
        ISSET2=1
        ENDIF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1-R1**2-R2**2)
        RFO4=128*(1+R1**2+R2**2-X1-X2)*(X1+X2)/(9*(-2+X1+X2)**2)-32*(-1
     &       +R1**4-6*R1**2*R2**2+R2**4+X1-R1**2*X1+3*R2**2*X1+X2
     &       +R1**2*X2-R2**2*X2-X1*X2)/(-1-R1**2+R2**2+X1)**2
     &       -32*(R1**2+R1**4+R2**2-6*R1**2*R2**2+R2**4-R1**2*X1
     &       +2*R2**2*X1+2*R1**2*X2-R2**2*X2-X1*X2)/
     &       ((-1-R1**2+R2**2+X1)*(-1+R1**2-R2**2+X2))-128*(-1+R1**4
     &       -6*R1**2*R2**2+R2**4+X1-R1**2*X1+R2**2*X1+X2+3*R1**2*X2
     &       -R2**2*X2-X1*X2)/(9*(-1+R1**2-R2**2+X2)**2)
     &       +16*(-1+R1**4-2*R2**2-R2**4-2*R1**2*X1+2*R2**2*X1+X1**2
     &       +X2-3*R1**2*X2+R2**2*X2+X1*X2)/
     &       ((-1-R1**2+R2**2+X1)*(-2+X1+ X2))
        RFO4=RFO4+16*(-1-2*R1**2-R1**4+R2**4+X1+R1**2*X1-3*R2**2*X1
     &       +2*R1**2*X2-2*R2**2*X2+X1*X2+X2**2)/
     &       (9*(1-R1**2+R2**2-X2)*(-2+X1+X2))
        RFO4=9D0*RFO4/128D0
        ISSET4=1
        ENDIF

C...q -> ~q ~g.
      ELSEIF(ICLASS.EQ.15) THEN
        IF(ICOMBI.EQ.1.OR.ICOMBI.EQ.3) THEN
        RLO1=PS*(1D0-R1**2+R2**2+2D0*R2)
        RFO1=32*(2*R2+X2)*(-1-R1**2-R2**2+X2)/(9*(-1+R1**2-R2**2+X2)**2)
     &       +8*(-1-R1**2-2*R1**2*R2-R2**2-2*R2**3+X1+R2*X1+R2**2*X1
     &       +3*X2/2-R1**2*X2/2+R2*X2-R2**2*X2/2-X1*X2/2)/
     &       ((-1-R1**2+R2**2+X1)*(-1+R1**2-R2**2+X2))+8*(2+2*R1**2-2*R2
     &       -2*R1**2*R2-6*R2**2-2*R2**3-3*X1-R1**2*X1+2*R2*X1
     &       +3*R2**2*X1+X1**2-X2-R1**2*X2+R2**2*X2+X1*X2)/
     &       (-1-R1**2+R2**2+X1)**2+32*(4+4*R1**2-4*R2**2-5*X1
     &       -R1**2*X1-2*R2*X1+R2**2*X1+X1**2-3*X2-R1**2*X2-2*R2*X2
     &       +R2**2*X2+X1*X2)/(9*(-2+X1+X2)**2)
        RFO1=RFO1+8*(3+3*R1**2-R2+R1**2*R2-5*R2**2-R2**3-4*X1-R1**2*X1
     &       +2*R2**2*X1+X1**2-2*X2-R2*X2+R2**2*X2+X1*X2)/
     &       ((-1-R1**2+R2**2+X1)*(2-X1-X2))+8*(-1-R1**2+R2+R1**2*R2
     &       -R2**2-R2**3+X1+R2*X1+R2**2*X1+2*X2+R1**2*X2-X1*X2/2
     &       -X2**2/2)/(9*(2-X1-X2)*(-1+R1**2-R2**2+X2))
        RFO1=9D0*RFO1/32D0
        ISSET1=1
        END IF
        IF(ICOMBI.EQ.2.OR.ICOMBI.EQ.3) THEN
        RLO2=PS*(1D0-R1**2+R2**2-2D0*R2)
        RFO2=32*(2*R2-X2)*(1+R1**2+R2**2-X2)/(9*(-1+R1**2-R2**2+X2)**2)
     &       +8*(-1-R1**2+2*R1**2*R2-R2**2+2*R2**3+X1-R2*X1+R2**2*X1
     &       +3*X2/2-R1**2*X2/2-R2*X2-R2**2*X2/2-X1*X2/2)/
     &       ((-1-R1**2+R2**2+X1)*(-1+R1**2-R2**2+X2))+8*(2+2*R1**2+2*R2
     &       +2*R1**2*R2-6*R2**2+2*R2**3-3*X1-R1**2*X1-2*R2*X1
     &       +3*R2**2*X1+X1**2-X2-R1**2*X2+R2**2*X2+X1*X2)/
     &       (-1-R1**2+R2**2+X1)**2+8*(3+3*R1**2+R2-R1**2*R2-5*R2**2
     &       +R2**3-4*X1-R1**2*X1+2*R2**2*X1+X1**2-2*X2+R2*X2+R2**2*X2
     &       +X1*X2)/((-1-R1**2+R2**2+X1)*(2-X1-X2))
        RFO2=RFO2+32*(4+4*R1**2-4*R2**2-5*X1-R1**2*X1+2*R2*X1+R2**2*X1
     &       +X1**2-3*X2-R1**2*X2+2*R2*X2+R2**2*X2+X1*X2)/
     &       (9*(-2+X1+X2)**2)+8*(-1-R1**2-R2-R1**2*R2-R2**2+R2**3+X1
     &       -R2*X1+R2**2*X1+2*X2+R1**2*X2-X1*X2/2-X2**2/2)/
     &       (9*(2-X1-X2)*(-1+R1**2-R2**2+X2))
        RFO2=9D0*RFO2/32D0
        ISSET2=1
        END IF
        IF(ICOMBI.EQ.4) THEN
        RLO4=PS*(1D0-R1**2+R2**2)
        RFO4=64*X2*(-1-R1**2-R2**2+X2)/(9*(-1+R1**2-R2**2+X2)**2)
     &       +16*(-1-R1**2-R2**2+X1+R2**2*X1+3*X2/2-R1**2*X2/2
     &       -R2**2*X2/2-X1*X2/2)/
     &       ((-1-R1**2+R2**2+X1)*(-1+R1**2-R2**2+X2))+16*(3+3*R1**2
     &       -5*R2**2-4*X1-R1**2*X1+2*R2**2*X1+X1**2-2*X2+R2**2*X2
     &       +X1*X2)/((-1-R1**2+R2**2+X1)*(2-X1-X2))
     &       +64*(4+4*R1**2-4*R2**2-5*X1-R1**2*X1+R2**2*X1+X1**2-3*X2
     &       -R1**2*X2+R2**2*X2+X1*X2)/(9*(-2+X1+X2)**2)
        RFO4=RFO4+16*(2+2*R1**2-6*R2**2-3*X1-R1**2*X1+3*R2**2*X1+X1**2
     &       -X2-R1**2*X2+R2**2*X2+X1*X2)/(-1-R1**2+R2**2+X1)**2
     &       +16*(-1-R1**2-R2**2+X1+R2**2*X1+2*X2+R1**2*X2-X1*X2/2
     &       -X2**2/2)/(9*(2-X1-X2)*(-1+R1**2-R2**2+X2))
        RFO4=9D0*RFO4/64D0
        ISSET4=1
        END IF

C...g -> ~g ~g. Use (9/4)*eikonal. May be changed in the future.
      ELSEIF(ICLASS.EQ.16) THEN
        RLO=PS
        IF(ICOMBI.EQ.0.OR.ICOMBI.EQ.1) THEN
          ANUM=0D0
        ELSEIF(ICOMBI.EQ.2) THEN
          ANUM=(2D0-X1-X2)**2
        ELSEIF(ICOMBI.EQ.3) THEN
          ANUM=ALPCOR*(2D0-X1-X2)**2
        ELSE
          ANUM=0.5D0*(2D0-X1-X2)**2
        ENDIF
        RFO=PS*2D0*((X1+X2-1D0+ANUM-R1**2-R2**2)/
     &       ((1D0+R1**2-R2**2-X1)*(1D0+R2**2-R1**2-X2))-
     &       R1**2/(1D0+R2**2-R1**2-X2)**2-
     &       R2**2/(1D0+R1**2-R2**2-X1)**2)
        RFO=9D0*RFO/4D0
        ICOMBI=0
      ENDIF

C...Find relevant LO and FO expression.
      IF(ICOMBI.EQ.0) THEN
      ELSEIF(ICOMBI.EQ.1.AND.ISSET1.EQ.1) THEN
        RLO=RLO1
        RFO=RFO1
      ELSEIF(ICOMBI.EQ.2.AND.ISSET2.EQ.1) THEN
        RLO=RLO2
        RFO=RFO2
      ELSEIF(ICOMBI.EQ.3.AND.ISSET1.EQ.1.AND.ISSET2.EQ.1) THEN
        RLO=ALPCOR*RLO1+(1D0-ALPCOR)*RLO2
        RFO=ALPCOR*RFO1+(1D0-ALPCOR)*RFO2
      ELSEIF(ISSET4.EQ.1) THEN
        RLO=RLO4
        RFO=RFO4
      ELSEIF(ICOMBI.EQ.4.AND.ISSET1.EQ.1.AND.ISSET2.EQ.1) THEN
        RLO=0.5D0*(RLO1+RLO2)
        RFO=0.5D0*(RFO1+RFO2)
      ELSEIF(ISSET1.EQ.1) THEN
        RLO=RLO1
        RFO=RFO1
      ELSE
        CALL PYERRM(16,'(PYMAEL:) not implemented ME code')
        RLO=1D0
        RFO=0D0
      ENDIF

C...Output.
      PYMAEL=RFO/RLO

      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
!> @file
!! Include fortran file for halting program in case
!! of omp activated
!! This file is useful in case of non-thread safe routine
!! that will raise a severe error if called inside a parallel region
!! this file should be included as !$ include to be considered only for OMP pragmas
!! @author
!!    Copyright (C) 2012-2013 BigDFT group
!!    This file is distributed under the terms of the
!!    GNU General Public License, see ~/COPYING file
!!    or http://www.gnu.org/copyleft/gpl.txt .
!!    For the list of contributors, see ~/AUTHORS

  logical :: in_omp
  logical, external :: omp_in_parallel,omp_get_nested
  in_omp=omp_in_parallel() .or. omp_get_nested()

  !!disable everything if we are into a OMP section
  !!timing routines are not thread-safe
  if(in_omp) then
     print *,'#SEVERE: critical f_lib routines cannot be called from within a parallel region! Exiting...'
     call f_err_severe_internal()
  end if

max_stars_repo_name: adegomme/aiida-bigdft-plugin


text:
module Mod_supf_elmope_4rdY
   use typre
   use Mod_Memor
   use Mod_Listen
   use Mod_Mesh
   use Mod_SUPFractionalStep
   use Mod_Element
   use Mod_NavierStokesElement
   use Mod_SUPF_Element
   use Mod_ConvectiveElement
   use Mod_TimeIntegrator
   use Mod_php_SetTimeIntegrator
   use Mod_sup_elmdir
   use Mod_nsm_elmdir
   use Mod_php_Elmdir
   use Mod_ThreeFieldElement
   use Mod_SupExacso
   use Mod_SupOperations

   implicit none

   class(SUPFractionalStepProblem), pointer :: a
   type(SupExacso) :: exacso
   real(rp), parameter :: zensi = 0.0_rp

   !Pointers
   type :: PPointer
      procedure(), NOPASS, pointer :: PostGaussElmats => NULL()
      procedure(), NOPASS, pointer :: ComputeAdvectionVelocity  => NULL()
      procedure(), NOPASS, pointer :: ResProRHS => NULL()
      procedure(), NOPASS, pointer :: ExtrapolationTerms => NULL()
      procedure(), NOPASS, pointer :: TauConstitutive  => NULL()
      procedure(), NOPASS, pointer :: ElmatUV => NULL()
      procedure(), NOPASS, pointer :: RHSUV => NULL()
      procedure(), NOPASS, pointer :: ConvectiveTerms => NULL()
      procedure(), NOPASS, pointer :: PreAssembly   => NULL()

   end type
   type(PPointer) :: ProcPointer

   !Hooks
   type :: PHook
      procedure(), NOPASS, pointer :: Initializations => NULL()
      procedure(), NOPASS, pointer :: Gathers => NULL()
      procedure(), NOPASS, pointer :: OnIeltyChange => NULL()
      procedure(), NOPASS, pointer :: PreGauss => NULL()
      procedure(), NOPASS, pointer :: Interpolates => NULL()
      procedure(), NOPASS, pointer :: InGauss => NULL()
      procedure(), NOPASS, pointer :: Testf => NULL()
      procedure(), NOPASS, pointer :: InGaussElmats => NULL()
      procedure(), NOPASS, pointer :: Finalizations => NULL()
      procedure(), NOPASS, pointer :: PhysicalProp  => NULL()

    end type
   type(PHook) :: ProcHook

   class(FiniteElement), pointer :: e => NULL()
   integer(ip) :: ielem,nelem

   real(rp), allocatable :: elmat(:,:,:,:)
   real(rp), allocatable :: elrhs(:,:)
   real(rp)              :: dvol

   type(TimeIntegratorDt1) :: Integrator
   integer(ip)           :: nsteps
   real(rp)              :: LHSdtinv

   integer(ip)           :: ielty0 = 0 !Previous element type

   real(rp), allocatable :: elvel(:,:,:)
   real(rp), allocatable :: elsig(:,:,:)
   real(rp), allocatable :: elpre(:,:),elext(:)
   real(rp), allocatable :: elmuv(:,:,:,:)
   real(rp), allocatable :: wrmat1(:,:)
   real(rp), allocatable :: AGradV(:), testf(:)
   real(rp), allocatable :: gpvel(:,:),gpadv(:),gpadv2(:)
   real(rp), allocatable :: gpsig(:,:)

   real(rp)             ::  val_aux

   real(rp)    :: chale(2),timom,tidiv,tisig,dvolt0,dvolt1,dvolt2,dvolt3,gprhs(3)
   integer(ip) :: nmean,auxtens
   real(rp)    :: acden,acvis,beta,auxvis,lambda,auxG
   real(rp)    :: reyno
   real(rp)    :: gpvno,gpvno2,gppre(1),gppre2(1),divvel,dummr,auxp1,auxp2,auxp3,auxp4,gppredif(1)
   !Residual Projections
   real(rp), allocatable :: elrep(:,:),gprep(:)
   real(rp), allocatable :: grpre(:,:),grvel(:,:),grvel2(:,:)
   !Split Oss
   real(rp), allocatable :: gpconv(:),elrepconv(:,:)
   real(rp), allocatable :: gplapl(:),elreplapl(:,:)
   real(rp)              :: gpdiv(1)
   real(rp), allocatable :: elrepdiv(:,:)

   integer(ip) :: idime,itime,igaus,itest,auxntens,auxndim,auxGrad,kdime,ldime,aux12
   integer(ip) :: currentbvess,bcstar

   !Level Set
   integer(ip)              :: inode,ipoin
   integer(ip)              :: elemStatus,ngauss_minus,ngauss_plus,ngaus_total,poinStatus
   real(rp), allocatable    :: weigp(:),xloc(:,:)

   !todo multy materials
   integer(ip) :: imat=1

   !#$COMPOSEPROCS 50
#include "COMPOSEPROCS_POINTERS_50.i90"

contains

#include "COMPOSEPROCS_SUBROUTINES_50.i90"

   !----------------------------------------------------------
   !SetPointers
   subroutine SetPointers
      use typre
      implicit none

      !External Procedures
      procedure() :: NULLSUB

      integer(ip) :: kfl_nonlinear,nelty

      call ResetProcedureComposition

      !-----------------------------------------------------------
      !Defaults

      !Pointers
      ProcPointer%PostGaussElmats => PostGaussElmats
      ProcPointer%ComputeAdvectionVelocity => ComputeAdvectionVelocity
      ProcPointer%ExtrapolationTerms  => sup_ExtraExtrapoaltion2d
      ProcPointer%ResProRHS       => ResProRHS2d
      ProcPointer%ElmatUV         => elmat2d
      ProcPointer%RHSUV           => rhs2d
      ProcPointer%ConvectiveTerms => elmatConvective
      !PreAssembly to use in free-surface and in enriched elements
      ProcPointer%PreAssembly => NULLSUB
      !Hooks
      ProcHook%Initializations => NULLSUB
      ProcHook%Gathers         => NULLSUB
      ProcHook%OnIeltyChange   => NULLSUB
      ProcHook%PreGauss        => NULLSUB
      ProcHook%Interpolates    => NULLSUB
      ProcHook%InGauss         => NULLSUB
      ProcHook%InGaussElmats   => NULLSUB
      ProcHook%Testf           => NULLSUB
      ProcHook%Finalizations   => NULLSUB
      ProcHook%PhysicalProp    => NULLSUB
      !Viscoelastic Part
      ProcPointer%TauConstitutive => TauElastic

      auxtens=(e%ndime-1)*(e%ndime-1)+2

      !MaterialProperties is called Always
      call ConcatenateProcedures(ProcHook%PhysicalProp,MaterialProperties)

      !Advection velocity
      if(e%ndime==3)then
         ProcPointer%ConvectiveTerms => elmatConvective3d
      endif

      !Advection velocity
      if (a%kfl_advec == 0) then
         ProcPointer%ComputeAdvectionVelocity => NULLSUB
         ProcPointer%ConvectiveTerms => NULLSUB
      endif

      !-----------------------------------------------------------
      !Non-linear elements
      call a%Mesh%IsNonLinearElement(ielem,kfl_nonlinear)
      if (kfl_nonlinear == 1) then
         call ConcatenateProcedures(ProcHook%InGauss,InGaussNonLinear)
      endif

      !-----------------------------------------------------------
      !Three-dimensional case
      if(e%ndime==3)then
         ProcPointer%ElmatUV         => elmat3d
         ProcPointer%RHSUV           => rhs3d
      endif

      !-----------------------------------------------------------
      !ResidualProjection
      !OSS
      if (a%kfl_repro == 1 .or. a%kfl_repro==0) then
         call runend('ASGS and OSS in Viscoelastic fractional dont work yet')
       endif

      !-----------------------------------------------------------

      !Split-Oss
      if(a%kfl_repro == 2 .or. a%kfl_repro == 3)then
         call ConcatenateProcedures(ProcHook%Initializations,AllocOssSplit)
         call ConcatenateProcedures(ProcHook%Gathers,GatherOssSplit)
         call ConcatenateProcedures(ProcHook%Interpolates,InterpolateOssSplit)
         call ConcatenateProcedures(ProcHook%InGaussElmats,InGaussElmatsOssSplit)
         call ConcatenateProcedures(ProcHook%Finalizations,DeallocOssSplit)
      end if

      !-----------------------------------------------------------
      !OSS constitutive
      if(a%kfl_repro == 2)then
         if(e%ndime == 3)then
            ProcPointer%ResProRHS       => ResProRHS3d
         end if
         call ConcatenateProcedures(ProcHook%Initializations,AllocRep)
         call ConcatenateProcedures(ProcHook%Gathers,GatherRep)
         call ConcatenateProcedures(ProcHook%Interpolates,InterpolateRep)
         call ConcatenateProcedures(ProcHook%InGaussElmats,InGaussElmatsRep)
         call ConcatenateProcedures(ProcHook%Finalizations,DeallocRep)
      end if

      !--------------------------------------------------------------------------------
      !If more than one element type, then pointers should be reset when ielty changes!
      call a%Mesh%GetNelty(nelty)
      if (nelty > 1) ProcHook%OnIeltyChange => OnIeltyChange

      !--------------------------------------------------------------------
      !Level set
      if(a%kfl_colev==1)then
         call ConcatenateProcedures(ProcHook%Initializations,AllocLevelsetTF)
         call ConcatenateProcedures(ProcHook%Finalizations,DeallocLevelSetTF)
         call PrependProcedure(ProcHook%PreGauss,CutelementsTF)
         call PrependProcedure(ProcHook%PhysicalProp,FluidPropertiesTF)
         if(a%kfl_fsurf==1)then
            if(a%kfl_fsurfLapla==1 .or. a%kfl_fsurfLapla==2)then
               call PrependProcedure(ProcPointer%PreAssembly,ElmatsToLapla)
            end if
            call ConcatenateProcedures(ProcPointer%PreAssembly,FreeSurfMatsTF)
         end if

      end if

   end subroutine

   !AdvectionVelocity
   subroutine ComputeAdvectionVelocity
      implicit none

      gpadv  = gpvel(:,1)
      gpadv2 = gpvel(:,2)
   end subroutine

   !----------------------------------------------------------
   !PostGauss Matrices
   subroutine PostGaussElmats
      implicit none

   end subroutine

   !-------------------------------------------------------------------
   !Physical Properties
   subroutine MaterialProperties
      implicit none
      !----------------------------------------------------------------
      !Incremental strategy (Continuation method in terms of lambda)
      acvis  = a%MatProp(imat)%LawViParam(1)
      acden  = a%MatProp(imat)%densi
      beta   = a%MatProp(imat)%LawViParam(2)

      !incremental scheme for Weissenberg number
      call a%IncrementalLambda(imat,a%MatProp(imat)%LawViParam(3), lambda)

      auxG   = a%MatProp(imat)%LawViParam(4)/((1.0_rp-beta) + 0.00001_rp)

   end subroutine

  !----------------------------------------------------------
   !NonLinear Elements
   subroutine InGaussNonLinear
      implicit none

      call e%elmder
      call e%elmhes

      !DvolsToZero
      dvolt0=0.0_rp
      dvolt1=0.0_rp
      dvolt2=0.0_rp
   end subroutine

   subroutine OnIeltyChange
      implicit none

      if (e%ielty /= ielty0) then
         call SetPointers
         ielty0 = e%ielty
      endif
   end subroutine

   !-------------------------------------------------------------------
   !FOR RESIDUAL PROJECTION
   subroutine AllocRep
      implicit none

      auxndim=a%ndofn
      call a%Memor%alloc(a%ndofn,e%mnode,elrep,'elrep','supf_elmope_1st')
      call a%Memor%alloc(a%ndofn,gprep,'gprep','supf_elmope_1st')

   end subroutine

   subroutine DeallocRep
      implicit none

      call a%Memor%dealloc(a%ndofn,e%mnode,elrep,'elrep','supf_elmope_1st')
      call a%Memor%dealloc(a%ndofn,gprep,'gprep','supf_elmope_1st')

   end subroutine

   subroutine GatherRep
      implicit none

      call e%gather(a%ndofn,elrep,a%repro)
   end subroutine

   subroutine InterpolateRep
      implicit none

      !Interpolate
      call e%interpg(a%ndofn,elrep,gprep)
   end subroutine

   subroutine InGaussElmatsRep
      implicit none

     call ProcPointer%ResProRHS

   end subroutine

   subroutine ResProRHS2d
      implicit none

      !Compute contributions to RHS : Block C
      call supf_elmrhu_oss(e,auxtens,tisig,dvol,gprep,beta,elrhs)

   end subroutine

   subroutine ResProRHS3d
      implicit none

      !Compute contributions to RHS : Block C
      call supf_elmrhu_oss3d(e,auxtens,tisig,dvol,gprep,beta,elrhs)

   end subroutine

   !------------------------------------------------------------------------------------------------------------
   !Split-OSS

   subroutine AllocOssSplit
      implicit none

      !Matrices alloc
      call a%Memor%alloc(e%ndime,gpconv, 'gpconv','supf_elmope_1st')
      call a%Memor%alloc(e%ndime,e%mnode,elrepconv, 'elrepconv','supf_elmope_1st')

      call a%Memor%alloc(1,e%mnode,elrepdiv, 'elrepdiv','supf_elmope_1st')

   end subroutine

   subroutine DeallocOssSplit
      implicit none

      !Matrices alloc
      call a%Memor%dealloc(e%ndime,gpconv, 'gpconv','supf_elmope_1st')
      call a%Memor%dealloc(e%ndime,e%mnode,elrepconv, 'elrepconv','supf_elmope_1st')

      call a%Memor%dealloc(1,e%mnode,elrepdiv, 'elrepdiv','supf_elmope_1st')

   end subroutine

   subroutine GatherOssSplit
      implicit none

      call e%gather(e%ndime,elrepconv,a%reproUGradU)
      call e%gather(1_ip,elrepdiv,a%reproDivU)

   end subroutine

   subroutine InterpolateOssSplit
      implicit none

      !Interpolate
      call e%interpg(e%ndime,elrepconv,gpconv)
      call e%interpg(1_ip,elrepdiv,gpdiv)

   end subroutine

   subroutine InGaussElmatsOssSplit
      implicit none

      call supf_rhu_soss_step1(e,dvol,acden,tidiv,gpdiv,elrhs)
!       call supf_rhu_soss_conv(e,dvol,acden,timom,AGradV,gpconv,elrhs)

   end subroutine

   !--------------------------------------------------------------------------------
   ! Stabilization Parameters
   subroutine TauElastic
      use Mod_sup_ComputeTidivVE
      use Mod_sup_ComputeTisig
      implicit none

      !Compute the stability parameters

      call ComputeTau(e,acden,acvis,gpvno,a%staco,chale,timom)
      call sup_ComputeTidivVE(e,timom,a%staco,chale,tidiv)
      call sup_ComputeTiSig(e,lambda,acvis,gpvno2,grvel,a%staco,chale,auxG,a%PTT_model,gpsig,tisig)

   end subroutine

   !-------------------------------------------------------------------
   !Physical Properties
   subroutine ViscosityLaw
      implicit none

   end subroutine

   !------------------------------------------------------------------
   !Laplacian elemental subroutines
   subroutine InGaussElmatsNonLinear
      implicit none

   end subroutine

   !-------------------------------------------------------------------------------
   !extrapolation terms

   subroutine sup_ExtraExtrapoaltion2d
   implicit none

      !Compute contributions to RHS : Block U
      call supf_elmrhu_corr(e,dvol,acden,LHSdtinv,gpvel(:,2),gppre(1),gpsig(:,1),elrhs)

      !Second order fractional step
      !Extra terms from the extrapolation initial values
      call supf_elmrhu_corr_extra(e,dvol,gppre2(1),gpsig(:,2),elrhs)

   end subroutine

   !matrices

   subroutine elmat2d

     !Constitutive stabilization term
     call supf_elmbuv_estS(e,dvol,acvis,acden,LHSdtinv,tidiv,tisig,beta,grvel,elmat)

   end subroutine

   subroutine elmat3d

     !Constitutive stabilization term
     call supf_elmbuv_estS3d(e,dvol,acvis,acden,LHSdtinv,tidiv,tisig,beta,grvel,elmat)

   end subroutine

   subroutine rhs2d

      !Compute contributions to RHS : Block U
      call supf_elmrhuY(e,acden,timom,tisig,beta,acvis,dvol,gpvel(:,1),gpvel(:,2),grvel,grvel2,AGradV,elrhs)

   end subroutine

   subroutine rhs3d

      !Compute contributions to RHS : Block U
      call supf_elmrhuY3d(e,acden,timom,tisig,beta,acvis,dvol,gpvel(:,1),gpvel(:,2),grvel,grvel2,AGradV,elrhs)

   end subroutine

   subroutine elmatConvective

      !Compute contributions to elemental matrix : Block U,V
      call supf_elmbuv1(e,dvol,acvis,beta,acden,LHSdtinv,aGradV,timom,gpvel,grvel,elmat)
      call supf_rhu_soss_conv(e,dvol,acden,timom,AGradV,gpconv,elrhs)
      !Compute contributions to RHS : Block U
      call supf_elmrhuYConv(e,acden,timom,tisig,beta,acvis,dvol,gpvel(:,1),gpvel(:,2),grvel,grvel2,AGradV,elrhs)
   end subroutine

   subroutine elmatConvective3d

      !Compute contributions to elemental matrix : Block U,V
      call supf_elmbuv13d(e,dvol,acvis,beta,acden,LHSdtinv,aGradV,timom,gpvel,grvel,elmat)
      call supf_rhu_soss_conv(e,dvol,acden,timom,AGradV,gpconv,elrhs)
      !Compute contributions to RHS : Block U
      call supf_elmrhuYConv3d(e,acden,timom,tisig,beta,acvis,dvol,gpvel(:,1),gpvel(:,2),grvel,grvel2,AGradV,elrhs)

   end subroutine

   !-------------------------------------------------------
   !LevelSet

   subroutine AllocLevelsetTF
      implicit none
      integer(ip) :: auxdim,auxdim1

      auxdim=(e%ndime+1)
      auxdim1=(e%ndime-1)*(e%ndime-1)+2

      call a%Memor%alloc(auxdim*auxdim1,weigp,'weigp','supf_elmope_1st')
      call a%Memor%alloc(e%ndime,auxdim*auxdim1,xloc,'xloc','supf_elmope_1st')

   end subroutine

   subroutine DeallocLevelSetTF
      implicit none
      integer(ip) :: auxdim,auxdim1

      auxdim=(e%ndime+1)
      auxdim1=(e%ndime-1)*(e%ndime-1)+2

      call a%Memor%dealloc(auxdim*auxdim1,weigp,'weigp','supf_elmope_1st')
      call a%Memor%dealloc(e%ndime,auxdim*auxdim1,xloc,'xloc','supf_elmope_1st')

   end subroutine

   subroutine CutelementsTF
      implicit none
      integer(ip) :: elemStatus

      call a%CutMesh%GetElementType(ielem,elemStatus)

      if(elemStatus ==0)then

         ngaus_total=0

         call a%CutMesh%GetNgaussSide(ielem,ngauss_plus,ngauss_minus)
         call a%CutMesh%GetWeigpCut(ielem,e%ndime,weigp)
         call a%CutMesh%GetXlocCut(ielem,e%ndime,xloc)

         ngaus_total = ngauss_minus+ngauss_plus

         !The rutine give the needed shape functions associated
         if(a%kfl_fsurf==1) weigp(1:ngauss_minus)=0.0_rp

         call e%SetParticularGaussPoints(a%Memor,ngaus_total,xloc,weigp(:))

      end if

   end subroutine

   subroutine FluidPropertiesTF
      implicit none
      integer(ip) :: elemStatus

      call a%CutMesh%GetElementType(ielem,elemStatus)

      if(elemStatus==1)then
         imat=1
         acden=a%MatProp(imat)%densi
         acvis=a%MatProp(imat)%visco
      elseif(elemStatus==-1)then
         imat=2
         acden=a%MatProp(imat)%densi
         acvis=a%MatProp(imat)%visco
      elseif(elemStatus ==0)then

         ngauss_minus=0
         ngauss_plus=0
         call a%CutMesh%GetNgaussSide(ielem,ngauss_plus,ngauss_minus)

         if(e%igaus<=ngauss_minus)then
            imat=2
            acden=a%MatProp(imat)%densi
            acvis=a%MatProp(imat)%visco
         elseif(e%igaus>ngauss_minus)then
            imat=1
            acden=a%MatProp(imat)%densi
            acvis=a%MatProp(imat)%visco
         end if

      end if

   end subroutine

   subroutine FreeSurfMatsTF
      implicit none
      integer(ip) :: elemStatus

      call a%CutMesh%GetElementType(ielem,elemStatus)

      if(elemStatus==-1)then
         do inode=1,e%pnode
            ipoin=e%lnods(inode)
            call a%CutMesh%GetPointType(ipoin,poinStatus)
            if(poinStatus==-1)then

               elmat(:,inode,:,1:e%pnode) = 0.0_rp
               elrhs(:,inode)=0.0_rp

            end if
         end do
      end if
   end subroutine

   subroutine ElmatsToLapla
      implicit none
      integer(ip)   :: elemStatus

      call a%CutMesh%GetElementType(ielem,elemStatus)
      if(elemStatus==-1)then

         elmat = 0.0_rp
         elrhs = 0.0_rp

         do inode=1,e%pnode
            do idime=1,e%ndime
               elmat(idime,inode,idime,inode) = elmat(idime,inode,idime,inode) + 1.0_rp

            end do
         end do

         elrhs(1:e%ndime,:) = elvel(1:e%ndime,:,1)

      end if

   end subroutine

end module

subroutine supf_elmope_4rdY(SUPFProblem)
   use Mod_SUPFractionalStep
   use Mod_supf_elmope_4rdY
   implicit none
   class(SUPFractionalStepProblem), target :: SUPFProblem

   a=>SUPFProblem

   call a%Mesh%ElementAlloc(e,a%Memor,'DefaultRule','supf_elmope_1st')

   !Set the Pointers for execution
   ielem = 1
   call SetPointers

   !Set Time Integrator
   call php_SetTimeIntegrator(a,Integrator,LHSDtinv,nsteps)

   !Matrices Alloc
   call a%Memor%alloc(e%ndime,e%mnode,e%ndime,e%mnode,elmat,'elmat','supf_elmope_1st')
   call a%Memor%alloc(e%ndime,e%mnode,elrhs,'elrhs','supf_elmope_1st')
   call a%Memor%alloc(e%mnode,e%mnode,wrmat1,'wrmat1','supf_elmope_1st')
   call a%Memor%alloc(e%ndime,e%mnode,e%ndime,e%mnode,elmuv,'elmuv','supf_elmope_1st')

   !Other arrays alloc
   call a%Memor%alloc(e%ndime,e%mnode,3,elvel,'elvel','supf_elmope_1st')
   call a%Memor%alloc(auxtens,e%mnode,2,elsig,'elsig','supf_elmope_1st')
   call a%Memor%alloc(      e%mnode,2,elpre,'elpre','supf_elmope_1st')
   call a%Memor%alloc(e%ndime,3,gpvel,'gpvel','supf_elmope_1st')
   call a%Memor%alloc(auxtens,2,gpsig,'gpsig','supf_elmope_1st')
   call a%Memor%alloc(e%ndime,gpadv,'gpadv','supf_elmope_1st')
   call a%Memor%alloc(e%ndime,gpadv2,'gpadv2','supf_elmope_1st')
   call a%Memor%alloc(e%ndime,elext,'elext','supf_elmope_1st')
   call a%Memor%alloc(e%mnode,AGradV,'AGradV','supf_elmope_1st')
   call a%Memor%alloc(e%mnode,testf,'testf','supf_elmope_1st')

   !gradients
   call a%Memor%alloc(e%ndime,e%ndime,grvel,'grvel','supf_elmope_1st')
   call a%Memor%alloc(e%ndime,e%ndime,grvel2,'grvel2','supf_elmope_1st')

   !Statistics
   call a%InitStats

   !Hook
   call ProcHook%Initializations

   call a%Mesh%GetNelem(nelem)
   elements : do ielem = 1,nelem
      !Load Element
      call a%Mesh%ElementLoad(ielem,e)

      !Hook
      call ProcHook%OnIeltyChange

      !Hook
      call ProcHook%PreGauss

      !Initializations
      elmat=0.0_rp
      elrhs=0.0_rp
      wrmat1 = 0.0_rp

      !Gathers
      !iteration i-1
      call e%gather(e%ndime,elvel(:,:,1),a%veloc(:,:,1))
      !intermediate values
      call e%gather(e%ndime,elvel(:,:,2),a%veloc(:,:,2))
      call e%gather(auxtens,elsig(:,:,1),a%sigma(:,:,2))
      call e%gather(1,elpre(:,1),a%press(:,1))
      !time step n
      call e%gather(e%ndime,elvel(:,:,3),a%veloc(:,:,3))
      call e%gather(auxtens,elsig(:,:,2),a%sigma(:,:,3))
      call e%gather(1,elpre(:,2),a%press(:,3))

      !Hook
      call ProcHook%Gathers

      !Cartesian derivatives and Jacobian at center of gravity
      call e%elmdcg

      !Element length at center of gravity
      call e%elmlen

      ! Compute the characteristic length chale
      call elmchl(e,a%kfl_advec,elvel,chale)

      !DvolsToZero
      dvolt0=0.0_rp
      dvolt1=0.0_rp
      dvolt2=0.0_rp

      !Gauss Point Loop
      gauss_points : do igaus=1,e%pgaus
         e%igaus = igaus

         !Hook
         call ProcHook%InGauss

         dvol = e%weigp(e%igaus)*e%detjm

         !Interpolate
         !iteration i-1
         call e%interpg(e%ndime,elvel(:,:,1),gpvel(:,1)) !i-1 iteration
         !intermediate values
         call e%interpg(e%ndime,elvel(:,:,2),gpvel(:,2)) !j-1 time
         call e%interpg(auxtens,elsig(:,:,1),gpsig(:,1)) !i-1 iteration
         call e%interpg(1,elpre(:,1),gppre(1))
         !time step n
         call e%interpg(e%ndime,elvel(:,:,3),gpvel(:,3)) !i-1 iteration
         call e%interpg(auxtens,elsig(:,:,2),gpsig(:,2)) !i-1 iteration
         call e%interpg(1,elpre(:,2),gppre2(1))

         !velocity and pressure gradient
         call e%gradient(e%ndime,elvel(:,:,1),grvel)

         !velocity and pressure gradient
         call e%gradient(e%ndime,elvel(:,:,2),grvel2)

        !Hook
         call ProcHook%Interpolates

         !Physical Parameters
         call a%GetPhysicalParameters(imat,acden,acvis)

         !Hook
         call ProcHook%PhysicalProp

         !Advection velocity
         call ProcPointer%ComputeAdvectionVelocity

         !Advection velocity norm
         call vecnor(gpadv,e%ndime,gpvno,2)
         !Advection velocity norm used in alfa3
         call vecnor(gpvel(:,1),e%ndime,gpvno2,2)

         !Compute a·grad(V)
         call ComputeAGradV(e,gpadv,AGradV)

         call ProcPointer%TauConstitutive

         elext =0.0_rp

         !Adjoint Test Function
         !Stabilization terms : -tau L*v
         call nsm_ComputeTestf(e,acden,timom,AGradV,testf)
         !Hook
         call ProcHook%Testf

         !Volumes Times Taus
         dvolt0=dvol       + dvolt0                ! w(gp)*detjm
         dvolt1=dvol*timom + dvolt1                ! w(gp)*detjm*tau1
         dvolt2=dvol*tidiv + dvolt2                ! w(gp)*detjm*tau2

         !----------------------------------------------------------------------
         ! Elemental matrix elmat and RHS

         call ProcPointer%ElmatUV
         call ProcPointer%ConvectiveTerms
         call ProcPointer%RHSUV

         !--------------------------------------------------------------------
         ! Second order fractional step
         !Terms involving the extrapolate pressure and stress to the RHS
         call ProcPointer%ExtrapolationTerms

         !Hook
         call ProcHook%InGaussElmats

         !Statistics
         call a%InGaussStats(acden,acvis,gpvno,chale,timom)

      enddo gauss_points

      !Post Gauss Matrices
      !Pointer
      call ProcPointer%PostGaussElmats

      !Matrix composition
      do idime = 1,e%ndime
         elmat(idime,1:e%pnode,idime,1:e%pnode) = elmat(idime,1:e%pnode,idime,1:e%pnode) + wrmat1(1:e%pnode,1:e%pnode)
      enddo

      !Pre Assembly Modifications
      !Pointer
      call ProcPointer%PreAssembly

      !ACA TENGO DUDAS CON LAS CONDICIONES DE BORDE
      !Dirichlet Boundary Conditions
      call nsm_rotdir(a,e,e%ndime,elmat,elrhs)
      !inside the subroutine
      ! php_elmdir(a,e,ndofn,ndofbc,ndofbcstart,currentbvess,elmat,elrhs)
      currentbvess=auxtens+1
      bcstar=0_ip

      if(a%kfl_confi==1.and.a%MatProp(imat)%lawvi<0)then
         call sup_Exaelmdir(a,e,e%ndime,e%ndime,bcstar,currentbvess,elmat,elrhs)
      elseif(a%kfl_confi==1.and.a%kfl_bc_number>0)then
         call sup_Exaelmdir(a,e,e%ndime,e%ndime,bcstar,currentbvess,elmat,elrhs)
      else
         call php_elmdir(a,e,e%ndime,e%ndime,bcstar,currentbvess,elmat,elrhs)
      end if

      !Assembly
      call a%LinearSystem%Assembly(e,elmat,elrhs)

   enddo elements

   call a%FinalizeStats

   !Hook
   call ProcHook%Finalizations

   !Matrices deAlloc
   call a%Memor%dealloc(e%ndime,e%mnode,e%ndime,e%mnode,elmat,'elmat','supf_elmope_1st')
   call a%Memor%dealloc(e%ndime,e%mnode,elrhs,'elrhs','supf_elmope_1st')
   call a%Memor%dealloc(e%mnode,e%mnode,wrmat1,'wrmat1','supf_elmope_1st')
   call a%Memor%dealloc(e%ndime,e%mnode,e%ndime,e%mnode,elmuv,'elmuv','supf_elmope_1st')

   !Other arrays dealloc
   call a%Memor%dealloc(e%ndime,e%mnode,3,elvel,'elvel','supf_elmope_1st')
   call a%Memor%dealloc(auxtens,e%mnode,2,elsig,'elsig','supf_elmope_1st')
   call a%Memor%dealloc(      e%mnode,2,elpre,'elpre','supf_elmope_1st')
   call a%Memor%dealloc(e%ndime,3,gpvel,'gpvel','supf_elmope_1st')
   call a%Memor%dealloc(auxtens,2,gpsig,'gpsig','supf_elmope_1st')
   call a%Memor%dealloc(e%ndime,gpadv,'gpadv','supf_elmope_1st')
   call a%Memor%dealloc(e%ndime,gpadv2,'gpadv2','supf_elmope_1st')
   call a%Memor%dealloc(e%ndime,elext,'elext','supf_elmope_1st')
   call a%Memor%dealloc(e%mnode,AGradV,'AGradV','supf_elmope_1st')
   call a%Memor%dealloc(e%mnode,testf,'testf','supf_elmope_1st')

   !gradients dealloc
   call a%Memor%dealloc(e%ndime,e%ndime,grvel,'grvel','supf_elmope_1st')
   call a%Memor%dealloc(e%ndime,e%ndime,grvel2,'grvel2','supf_elmope_1st')

   !ElementDealloc
   call a%Mesh%ElementDealloc(e,a%Memor,'DefaultRule','supf_elmope_1st')

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM