text:
KPFA was founded in 1949 and is located in Berkeley, but simulcasts on two other stations, one in Berkeley and one is Fresno. Combined, their broadcast reaches 1/3 of the state, according to their website. Their programming consists of eclectic music, public affairs and news coverage. It is a Pacifica station with a mission:

     To promote cultural diversity and pluralistic community expression
     To contribute to a lasting understanding between individuals of all nations, races, creeds and colors
     To promote freedom of the press and serve as a forum for various viewpoints
     To maintain an independent funding base

KDVS likes to consider itself closer in kin to KPFA than to UC Berkeleys campus station, http://kalx.berkeley.edu KALX.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    @(#)chk_tiedta.f	20.3 12/21/96
      integer function chk_tiedta(ptr, field, count, out_buffer)
      integer ptr, field, count
      character out_buffer(10)*120

C     This subroutine checks ARCINP(*, PTR) extensively for data errors.

      include 'ipfinc/parametr.inc'
      include 'ipfinc/arcntl.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/prt.inc'

      chk_tiedta = 0
      if ((field .eq. 0 .or. field .eq. 27) .and.
     &    (arcinp(ptr) .lt. -6000.0 .or. arcinp(ptr) .gt. 6000.0)) then
        chk_tiedta = 1
        count = min0 (count + 1, 10)
        write (out_buffer(count), 10020) arcinp(ptr)
10020   format (' Export (', f10.1, ') < -6000 or > 6000 ')
      endif

  220 return
      end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      integer function sysruk (task, cmd, rukarf, rukint)
      integer rukarf
      integer rukint
      integer*2 task(*)
      integer*2 cmd(*)
      integer i
      integer ntasks
      integer lmarg
      integer rmarg
      integer maxch
      integer ncol
      integer rukean
      integer envgei
      integer envscn
      logical streq
      logical xerpop
      logical xerflg
      common /xercom/ xerflg
      integer iyy
      integer dp(2)
      integer*2 dict(7)
      integer*2 st0001(9)
      integer*2 st0002(6)
      integer*2 st0003(3)
      integer*2 st0004(6)
      integer*2 st0005(6)
      integer*2 st0006(4)
      integer*2 st0007(6)
      integer*2 st0008(2)
      integer*2 st0009(29)
      integer*2 st0010(25)
      save
      data dict /119, 99,115,112,105,120, 0/
      data (st0001(iyy),iyy= 1, 8) /116,116,121,110, 99,111,108,115/
      data (st0001(iyy),iyy= 9, 9) / 0/
      data st0002 / 99,104,100,105,114, 0/
      data st0003 / 99,100, 0/
      data st0004 /104,111,109,101, 36, 0/
      data st0005 / 72, 79, 77, 69, 36, 0/
      data st0006 /115,101,116, 0/
      data st0007 /114,101,115,101,116, 0/
      data st0008 / 9, 0/
      data (st0009(iyy),iyy= 1, 8) /105,110,118, 97,108,105,100, 32/
      data (st0009(iyy),iyy= 9,16) /115,101,116, 32,115,116, 97,116/
      data (st0009(iyy),iyy=17,24) /101,109,101,110,116, 58, 32, 39/
      data (st0009(iyy),iyy=25,29) / 37,115, 39, 10, 0/
      data (st0010(iyy),iyy= 1, 8) /105,110,118, 97,108,105,100, 32/
      data (st0010(iyy),iyy= 9,16) / 83, 69, 84, 32,105,110, 32, 73/
      data (st0010(iyy),iyy=17,24) / 82, 65, 70, 32, 77, 97,105,110/
      data (st0010(iyy),iyy=25,25) / 0/
      data (dp(iyy),iyy= 1, 2) / 1, 0/
      data lmarg /5/, maxch /0/, ncol /0/, rukean /3/
      data ntasks /0/
         if (.not.(ntasks .eq. 0)) goto 110
            i=1
120         if (.not.(dp(i) .ne. 0)) goto 122
121            i=i+1
               goto 120
122         continue
            ntasks = i - 1
110      continue
         if (.not.(task(1) .eq. 63)) goto 130
            call xerpsh
            rmarg = envgei (st0001)
            if (.not.xerpop()) goto 140
               rmarg = 80
140         continue
            call strtbl (4, dict, dp, ntasks, lmarg, rmarg, maxch, ncol)
            sysruk = (0)
            goto 100
130      continue
         if (.not.(streq(task,st0002) .or. streq(task,st0003))) goto 150
            call xerpsh
            if (.not.(cmd(rukarf) .eq. 0)) goto 170
               call xerpsh
               call xfchdr(st0004)
               if (.not.xerpop()) goto 180
                  call xfchdr(st0005)
180            continue
               goto 171
170         continue
               call xfchdr(cmd(rukarf))
171         continue
162         if (.not.xerpop()) goto 160
               if (.not.(rukint .eq. 1)) goto 190
                  call erract (rukean)
                  if (xerflg) goto 100
                  goto 191
190            continue
191            continue
160         continue
            sysruk = (0)
            goto 100
150      continue
         if (.not.(streq(task,st0006) .or. streq(task,st0007))) goto 200
            call xerpsh
            if (.not.(cmd(rukarf) .eq. 0)) goto 220
               call envlit (4, st0008, 1)
               call xffluh(4)
               goto 221
220         continue
            if (.not.(envscn (cmd) .le. 0)) goto 230
               if (.not.(rukint .eq. 1)) goto 240
                  call eprinf (st0009)
                  call pargsr (cmd)
                  goto 241
240            continue
                  goto 91
241            continue
230         continue
221         continue
212         if (.not.xerpop()) goto 210
               if (.not.(rukint .eq. 1)) goto 250
                  call erract (rukean)
                  if (xerflg) goto 100
                  goto 251
250            continue
91                call syspac (0, st0010)
251            continue
210         continue
            sysruk = (0)
            goto 100
200      continue
151      continue
131      continue
         if (.not.(streq (task, dict(dp(1))))) goto 260
            call twcspx
            sysruk = (0)
            goto 100
260      continue
         sysruk = (-1)
         goto 100
100      return
      end
c     rukint  ruk_interact
c     sysruk  sys_runtask
c     envscn  envscan
c     twcspx  t_wcspix
c     envgei  envgeti
c     syspac  sys_panic
c     eprinf  eprintf
c     rukarf  ruk_argoff
c     rukean  ruk_eawarn
c     pargsr  pargstr
c     envlit  envlist

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: ysBach/irafdocgen

max_forks_repo_name: ysBach/irafdocgen


text:
module type_tutorial_step
#include <messenger.h>
    implicit none

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   5/17/2017
    !!
    !!
    !-----------------------------------------------------------------------------------
    type, public, abstract :: tutorial_step_t

        character(:),               allocatable :: title
        !class(tutorial_action_t),   allocatable :: action

    contains

        procedure           :: execute
        !procedure(deferred) :: execute

        procedure           :: set_title
        procedure           :: get_title

    end type tutorial_step_t
    !***********************************************************************************

contains

    !>  Set the title of the step.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   5/17/2017
    !!
    !-----------------------------------------------------------------------------------
    subroutine set_title(self,title)
        class(tutorial_step_t), intent(inout)   :: self
        character(*),           intent(in)      :: title

        self%title = title

    end subroutine set_title
    !***********************************************************************************

    !>  Return the title of the step.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   5/17/2017
    !!
    !-----------------------------------------------------------------------------------
    function get_title(self) result(title)
        class(tutorial_step_t), intent(in)  :: self

        character(:),   allocatable :: title

        title = self%title

    end function get_title
    !***********************************************************************************

    !>  Execute the step procedures.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   5/17/2017
    !!
    !-----------------------------------------------------------------------------------
    subroutine execute(self)
        class(tutorial_step_t), intent(inout)   :: self

    end subroutine execute
    !***********************************************************************************

end module type_tutorial_step

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
program ch0811
    ! Rank 2 Arrays and the Sum Intrinsic Funcoltion
    implicit none

    integer, parameter :: nrow = 5
    integer, parameter :: ncol = 6
    real, dimension(1:nrow*ncol) :: results = [ &
        50., 47., 28., 89., 30., 46., 37., 67., 34., 65., 68., &
        98., 25., 45., 26., 48., 10., 36., 89., 56., 33., 45., &
        30., 65., 68., 78., 38., 76., 98., 65. ]
    real, dimension(1:nrow,1:ncol) :: exam_results = 0.
    real, dimension(1:nrow) :: people_average = 0.
    real, dimension(1:ncol) :: subject_average = 0.

    exam_results = reshape(results, [nrow,ncol], [0.,0.], [2,1])
    exam_results(1:nrow,3) = exam_results(1:nrow,3)*2.5
    people_average = sum(exam_results, 2)
    people_average = people_average/ncol
    subject_average = sum(exam_results, 1)
    subject_average = subject_average/nrow
    print *, ' People averages'
    print *, people_average
    print *, ' Subject averages'
    print *, subject_average
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
Judy Sakaki was the UCD Administrative Officers Vice Chancellor of http://studentaffairs.ucdavis.edu Student Affairs from July 2002 until January 2007.  This entailed a lot of bureaucratic stuff in Mrak Hall, but shell still probably payed more attention to you than, say, Stan Nosek, whos in charge of administration.  Her office was in 476 Mrak Hall.  Shes got a little http://studentleadership.ucdavis.edu/focus_campusleaders/judy_sakaki.html profile c/o some sort of leadership website.

Sakaki indirectly oversaw many aspects of student life, including Campus Unions, UCD Athletics Athletics, Advising Services, the UCD Bookstore, Financial Aid, the Registrar, Student Housing, Cowell Student Health Center, and Undergraduate Admissions.

Sakaki was also the administration goto person on issues of Cow Tipping cow tipping.

Judy used to send you a card before your 21st birthday telling you not to die because it makes the UC Davis University look bad.

She is now the Vice President for Student Affairs at the UC Office of the President.

She didnt send me a card. :(  Which is too bad, since that was a point in my life where I really wanted to die. Users/BrentLaabs

  Im guessing you turned 21 before Tragic Events David Thorton.  Users/ArlenAbraham


text:
Name: John Terning

Phone Number: (530) 7522749

Office: 435 Physics and Geology Physics, UC Davis

Im an Associate Professor in Physics; I arrived at UCD in Jan. 2005.
My http://particle.physics.ucdavis.edu/john/index.php research concentrates on particle physics and the origin of mass.


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

! ! in this section you can have comments (which need a second !) or
! ! lines with two words on them that begins QTY_xxx "units"
! ! can have trailing comments after:  QTY_xxx  "units" ! comment
!
! BEGIN DART PREPROCESS QUANTITY DEFINITIONS
!
!
! ! kinds for CICE added by C. Bitz
!   QTY_SEAICE_AGREG_CONCENTR           "none"
!   QTY_SEAICE_AGREG_VOLUME             "none"
!   QTY_SEAICE_AGREG_SNOWVOLUME         "none"
!   QTY_SEAICE_AGREG_THICKNESS          "none"
!   QTY_SEAICE_AGREG_SNOWDEPTH          "none"
!   QTY_U_SEAICE_COMPONENT              "none"
!   QTY_V_SEAICE_COMPONENT              "none"
!   QTY_SEAICE_ALBEDODIRVIZ             "none"
!   QTY_SEAICE_ALBEDODIRNIR             "none"
!   QTY_SEAICE_ALBEDOINDVIZ             "none"
!   QTY_SEAICE_ALBEDOINDNIR             "none"
!   QTY_SEAICE_CONCENTR                 "none"
!   QTY_SEAICE_VOLUME                   "none"
!   QTY_SEAICE_SNOWVOLUME               "none"
!   QTY_SEAICE_SURFACETEMP              "none"
!   QTY_SEAICE_FIRSTYEARAREA            "none"
!   QTY_SEAICE_ICEAGE                   "none"
!   QTY_SEAICE_LEVELAREA                "none"
!   QTY_SEAICE_LEVELVOLUME              "none"
!   QTY_SEAICE_MELTPONDAREA             "none"
!   QTY_SEAICE_MELTPONDDEPTH            "none"
!   QTY_SEAICE_MELTPONDLID              "none"
!   QTY_SEAICE_MELTPONDSNOW             "none"
!   QTY_SEAICE_SALINITY001              "none"
!   QTY_SEAICE_SALINITY002              "none"
!   QTY_SEAICE_SALINITY003              "none"
!   QTY_SEAICE_SALINITY004              "none"
!   QTY_SEAICE_SALINITY005              "none"
!   QTY_SEAICE_SALINITY006              "none"
!   QTY_SEAICE_SALINITY007              "none"
!   QTY_SEAICE_SALINITY008              "none"
!   QTY_SEAICE_ICEENTHALPY001           "none"
!   QTY_SEAICE_ICEENTHALPY002           "none"
!   QTY_SEAICE_ICEENTHALPY003           "none"
!   QTY_SEAICE_ICEENTHALPY004           "none"
!   QTY_SEAICE_ICEENTHALPY005           "none"
!   QTY_SEAICE_ICEENTHALPY006           "none"
!   QTY_SEAICE_ICEENTHALPY007           "none"
!   QTY_SEAICE_ICEENTHALPY008           "none"
!   QTY_SEAICE_SNOWENTHALPY001          "none"
!   QTY_SEAICE_SNOWENTHALPY002          "none"
!   QTY_SEAICE_SNOWENTHALPY003          "none"
!   QTY_SOM_TEMPERATURE                 "none"
!   QTY_SEAICE_FY                       "none"
!   QTY_SEAICE_AGREG_FY                 "none"
!   QTY_SEAICE_AGREG_SURFACETEMP        "none"
!   QTY_SEAICE_AGREG_FREEBOARD          "none"
!   QTY_SEAICE_CATEGORY                 "none"
!
! END DART PREPROCESS QUANTITY DEFINITIONS

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
      program drdilu
c>> 2001-05-22 DRDILU Krogh Minor change for making .f90 version.
c>> 1994-12-21 DRDILU  Krogh  Latest version.
c--D replaces "?": DR?ILU, ?ILUP
c Demonstration driver for DILUP.
c Given table of sin(x), x = 0, .5, 1., ..., interpolates for x =
c -.1, .1, and, 8.3 using polynomial degrees 2, 3, 4, ..., 10.  Errors
c are estimated in all cases.
c
      integer NTAB, NX, NDEG1, NDEG2
      parameter (NTAB=40, NX=3, NDEG1=2, NDEG2=10)
      double precision XT(2), YT(NTAB), X(NX), Y, EOPT(1), H, ANS
      parameter (H = .5D0)
      integer I, NDEG, IOPT(3), LUP
      data X / -.1D0, .1D0, 8.3D0 /
c                       Set IOPT to get an error estimate.
      data IOPT / 0, 1, 0 /
      data LUP / 3 /
c
c Compute the XT and YT tables
      XT(1) = 0.D0
      XT(2) = H
      do 10 I = 1, NTAB
         YT(I) = sin(dble(I-1) * H)
   10 continue
      print *, 'IOP(1) NDEG    X         Y     Est. Error True Error'
      do 30 I = 1, NX
         ANS = sin(X(I))
         do 20 NDEG = NDEG1, NDEG2
            call DILUP (X(I), Y, NTAB, XT, YT, NDEG, LUP, IOPT, EOPT)
            print '(I5, I5, F9.4, F12.8, 1P,E10.2, E11.2)',
     1       IOPT(1), NDEG, X(I), Y, EOPT(1), Y - ANS
   20    continue
      print *
   30 continue
      stop
      end


text:
!!
!! PARENT COMPONENT NAME NOT THE RENAMED PARENT TYPE
!!
!! The NAG compiler gives this error for the following example:
!!
!! $ nagfor nag-20180203c.f90
!! NAG Fortran Compiler Release 6.2(Chiyoda) Build 6205
!! Error: nag-bug-20180203c.f90, line 12: POINT is not a component of P
!!        detected at %@POINT
!!
!! I don't think this is correct. The name of the parent component is POINT
!! and not PUNKT. From the F08 standard:
!!
!! Note 4.51: "The name of the parent type might be a local name introduced
!! via renaming in a USE statement."  Here that is POINT.
!!
!! 4.5.7.2: "An extended type has a scalar, nonpointer, nonallocatable, parent
!! component with [...] The name of this component is the parent type name.
!!

module punkt_type
  type punkt
    real x, y
  end type
end module

use punkt_type, only: point => punkt
type, extends(point) :: color_point ! Note 4.51 (F08) POINT is the name of the parent type
  integer color
end type
type(color_point) :: p
p%point%x = 1.0   ! THE PARENT COMPONENT SHOULD BE P%POINT, NOT P%PUNKT
end

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20081120 17:57:48 nbsp Welcome to the Wiki Howdy Jef and Welcome to the Wiki!   My names Evan, pleased to meet you!  Thanks for adding a review to Round Table Pizza, but please dont delete the reviews of other editors, even the ones you disagree with.  I restored the reviews and photos.  Was there a reason you removed the photos as well?  Users/JabberWokky Evan JabberWokky Edwards


text:
C    @(#)opieqiv.f	20.4 2/13/96
      subroutine opieqiv (p, y, error)
      complex y(2,2)
      integer p, error
c
c     This subroutine computes the 2-port admittance matrix y
c     for obrnch(1,aq) from the ALTERNATE BASE CASE IN RESIDENCE!
c
      include 'ipfinc/parametr.inc'

      include 'ipfinc/alt_case.inc'
      include 'ipfinc/prt.inc'

c     set up record type code parameters...
      include 'ipfinc/brtype.inc'

      common /is_batch / is_batch

      complex cx
      integer q, aq
      character xbuf * 120
c
      error = 0
      k1 = okx(p)
      k2 = oky(p)
      ltype = obrtype(p)
      q = obrnch_ptr(p)
      aq = iabs(q)

      if (ltype.eq.BRTYP_L .or. ltype.eq.BRTYP_T .or.
     &    ltype.eq.BRTYP_TP .or. ltype.eq.BRTYP_E) then
c
c        compute transfer impedance
c
         r=obrnch(5,aq)
         x=obrnch(6,aq)
         z=r**2 + x**2
         if (z .eq. 0.0) then
            write (errbuf(1),110)
  110       format(' Branch has zero transfer impedance.')
            errbuf(2) = ' '
            call bcdbrn(p,xbuf)
            write (errbuf(3), 120) xbuf(1:80)
  120       format(13x,'(',a80,')')
            if (is_batch .eq. 0) then
               call prterx ('E',3)
            else
               call prterx ('F',3)
            endif
            error = 1
            z=1.0
         endif

         if (ltype .eq. BRTYP_L) then
c
c           "L" branch
c
            y(1,2)=-cmplx (r/z, -x/z)
            y(2,1) = y(1,2)
            y(1,1)=cmplx (obrnch(7,aq),obrnch(8,aq))-y(1,2)
            y(2,2)=y(1,1)

         else if (ltype .eq. BRTYP_T) then
c
c           "T" branch
c
            if (q .gt. 0) then
               tk=obrnch(9,aq)/oldbase(k1)
               tm=obrnch(10,aq)/oldbase(k2)
            else
               tk=obrnch(10,aq)/oldbase(k1)
               tm=obrnch(9,aq)/oldbase(k2)
            endif
            if (amax1(tk,tm).ge.1.40 .or. amin1(tk,tm).le.0.60) then
               if (q .gt. 0) then
                  tk=obrnch(10,aq)/oldbase(k1)
                  tm=obrnch(9,aq)/oldbase(k2)
               else
                  tk=obrnch(9,aq)/oldbase(k1)
                  tm=obrnch(10,aq)/oldbase(k2)
               endif
               if (amax1(tk,tm).lt.1.40 .and. amin1(tk,tm).gt.0.60)
     &            then
                  write(errbuf(1),130)
  130             format(' Backwards transformer taps are reversed.')
                  errbuf(2) = ' '
                  call bcdbrn(p,xbuf)
                  write (errbuf(3),120) xbuf(1:80)
                  temp = obrnch(9,aq)
                  obrnch(9,aq) = obrnch(10,aq)
                  obrnch(10,aq) = temp
                  call bcdbrn(p,xbuf)
                  write (errbuf(4),120) xbuf(1:80)
                  call prterx ('W',4)
               else
                  write (errbuf(1),140)
  140             format(' Transformer taps are beyond 40% of base kv.')
                  errbuf(2) = ' '
                  call bcdbrn(p,xbuf)
                  write (errbuf(3),120) xbuf(1:80)
                  if (is_batch .eq. 0) then
                     call prterx ('E',3)
                  else
                     call prterx ('F',3)
                  endif
                  if (tk.eq.0.0) tk=1.0
                  if (tm.eq.0.0) tm=1.0
               endif
            endif
            tkm = 1.0 / (tk*tm)
            y(1,1) = cmplx (0.5*tkm*obrnch(7,aq),
     &                     -0.5*tkm*abs(obrnch(8,aq)))
            y(2,2) = y(1,1)
            y(1,2) = -cmplx(r/(z*tk*tm),-x/(z*tk*tm))
            y(2,1) = y(1,2)
            y(1,1) = y(1,1) + cmplx (r/(z*tk*tk),-x/(z*tk*tk))
            y(2,2) = y(2,2) + cmplx (r/(z*tm*tm), -x/(z*tm*tm))

         else if (ltype .eq. BRTYP_TP) then
c
c           "TP" branch
c
            if (q .gt. 0) then
               tk=obrnch(9,aq)
               tm=obrnch(10,aq)/oldbase(k2)
            else
               tk=-obrnch(9,aq)
               tm=obrnch(10,aq)/oldbase(k1)
            endif
            if (tm .eq. 0.0) tm=1.0
            if (tm .ge. 1.40 .or. tm .le. 0.60) then
               write (errbuf(1),150)
  150          format(' Phase shifter tap2 is beyond 40% of base kv.')
               errbuf(2) = ' '
               call bcdbrn(p,xbuf)
               write (errbuf(3),120) xbuf(1:80)
               call prterx ('W',3)
               tm=1.0
            endif
            if (tk .ge. 120.1 .or. tk .le. -120.1) then
               write (errbuf(1), 160)
  160          format(' Phase shifter has excessive angle ',
     &                '(> 120 degrees)')
               errbuf(2) = ' '
               call bcdbrn(p,xbuf)
               write (errbuf(3),120) xbuf(1:80)
               call prterx ('W',3)
            endif
            angle=0.0174532925*tk
            y(1,2)=cmplx(0.5*obrnch(7,aq),-0.5*abs(obrnch(8,aq)))
            y(1,1)=y(1,2)+cmplx(r/z,-x/z)
            y(2,2)=y(1,2)+cmplx(r/z,-x/z)/cmplx(tm**2,0.0)
            y(1,2)=-cmplx(cos(angle),
     &                    sin(angle))*cmplx(r/(z*tm),-x/(z*tm))
            y(2,1)=-cmplx(cos(angle),
     &                   -sin(angle))*cmplx(r/(z*tm),-x/(z*tm))
            if (q .lt. 0) then
               cx=y(2,2)
               y(2,2)=y(1,1)
               y(1,1)=cx
            endif

         else
c
c           "E" branch
c
            y(1,2)= -cmplx (r/z, -x/z)
            y(2,1)=y(1,2)
            if (q .gt. 0) then
               y(1,1)=cmplx(obrnch(7,aq),obrnch(8,aq))-y(1,2)
               y(2,2)=cmplx(obrnch(9,aq),obrnch(10,aq))-y(2,1)
            else
               y(2,2)=cmplx(obrnch(7,aq),obrnch(8,aq))-y(2,1)
               y(1,1)=cmplx(obrnch(9,aq),obrnch(10,aq))-y(1,2)
            endif

         endif

      else if (ltype.eq.BRTYP_PEQ) then
c
c        "pi-equivalent" branch
c
         if (q .gt. 0) then
            y(1,1) = cmplx(obrnch(4,aq),obrnch(5,aq))
            y(1,2) = cmplx(obrnch(6,aq),obrnch(7,aq))
            y(2,1) = cmplx(obrnch(8,aq),obrnch(9,aq))
            y(2,2) = cmplx(obrnch(10,aq),obrnch(11,aq))
         else
            y(2,2) = cmplx(obrnch(4,aq),obrnch(5,aq))
            y(2,1) = cmplx(obrnch(6,aq),obrnch(7,aq))
            y(1,2) = cmplx(obrnch(8,aq),obrnch(9,aq))
            y(1,1) = cmplx(obrnch(10,aq),obrnch(11,aq))
         endif

      else
c
c        error-illegal branch type calling "pieqiv"
c
         write (errbuf(1), 170)
  170    format(' Illegal branch type processed by "PIEQIV"')
         errbuf(2) = ' '
         call bcdbrn(p,xbuf)
         write (errbuf(3),120) xbuf(1:80)
         call prterx ('W',3)
         error = 1
      endif

      return
      end


text:
C--------------------------------------------------------------------------
        subroutine ftgtcl(iunit,colnum,datcod,repeat,width,status)

C       get the datatype of the column, as well as the vector
C       repeat count and (if it is an ASCII character column) the
C       width of a unit string within the column.  This supports the
C       TFORMn = 'rAw' syntax for specifying arrays of substrings.

C       iunit   i  Fortran i/o unit number
C       colnum  i  number of the column (first column = 1)

C       datcod  i  returned datatype code
C       repeat  i  number of elements in the vector column
C       width   i  width of unit string in character columns
C       status  i  returned error status
C
C       written by Wm Pence, HEASARC/GSFC, November 1994

        integer iunit,colnum,datcod,repeat,width,status

C       COMMON BLOCK DEFINITIONS:--------------------------------------------
        integer nb,ne
        parameter (nb = 20)
        parameter (ne = 200)
        integer bufnum,chdu,hdutyp,maxhdu,hdstrt,hdend,nxthdr,dtstrt
        integer nxtfld
        logical wrmode
        common/ft0001/bufnum(199),chdu(nb),hdutyp(nb),maxhdu(nb),
     &  wrmode(nb),hdstrt(nb,ne),hdend(nb),nxthdr(nb),dtstrt(nb),nxtfld
C       END OF COMMON BLOCK DEFINITIONS------------------------------------

        integer ibuff,dummy
        character keywrd*8,tform*24,comm*20

        if (status .gt. 0)return

C       construct the keyword name
        call ftkeyn('TFORM',colnum,keywrd,status)

C       get the keyword value
        call ftgkys(iunit,keywrd,tform,comm,status)
        if (status .gt. 0)then
            call ftpmsg('Could not read the '//keywrd//' keyword.')
            return
        end if

C       parse the keyword value
        ibuff=bufnum(iunit)
        if (hdutyp(ibuff) .eq. 1)then
C           this is an ASCII table
            repeat=1
            call ftasfm(tform,datcod,width,dummy,status)

        else if (hdutyp(ibuff) .eq. 2)then
C           this is a binary table
            call ftbnfm(tform,datcod,repeat,width,status)

        else
C           error: this HDU is not a table
            status=235
            return
        end if
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
  3rd Year Animal Science/Economics Major
  Currently working at Quickly
  Part of UC Davis HKSA student organization (eventhough im not cantonese)
  Come to our events please~


text:
Address(Zamora Lane) is a residential street in East Davis.

  San Rafael Street
  Mallorca Lane
  Santiago Street


!
!
!                             AMG4PSBLAS version 1.0
!    Algebraic Multigrid Package
!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
!    (C) Copyright 2021
!
!        Salvatore Filippone
!        Pasqua D'Ambra
!        Fabio Durastante
!
!    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.
!      3. The name of the AMG4PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    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 AMG4PSBLAS GROUP OR ITS 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.
!
!
!
!  The aggregator object hosts the aggregation method for building
!  the multilevel hierarchy.
!
module amg_s_base_aggregator_mod

  use amg_base_prec_type, only : amg_sml_parms, amg_saggr_data
  use psb_base_mod, only : psb_sspmat_type, psb_lsspmat_type, psb_s_vect_type, &
       & psb_s_base_vect_type, psb_slinmap_type, psb_spk_, &
       & psb_ls_csr_sparse_mat,  psb_ls_coo_sparse_mat, &
       & psb_s_csr_sparse_mat,  psb_s_coo_sparse_mat, &
       & psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
       & psb_erractionsave, psb_error_handler, psb_success_, psb_toupper
  !
  !
  !
  !> \class amg_s_base_aggregator_type
  !!
  !!  It is the data type containing the basic interface definition for
  !!  building a multigrid hierarchy by aggregation. The base object has no attributes,
  !!  it is intended to be essentially an abstract type.
  !!
  !!
  !!  type amg_s_base_aggregator_type
  !!  end type
  !!
  !!
  !!  Methods:
  !!
  !!  bld_tprol   -   Build a tentative prolongator
  !!
  !!  mat_bld     -   Build prolongator/restrictor and coarse matrix ac
  !!
  !!  mat_asb     -   Convert  prolongator/restrictor/coarse matrix
  !!                  and fix their descriptor(s)
  !!
  !!  update_next -   Transfer information to the next level; default is
  !!                  to do nothing, i.e. aggregators at different
  !!                  levels are independent.
  !!
  !!  default     -   Apply defaults
  !!  set_aggr_type - For aggregator that have internal options.
  !!  fmt         - Return a short string description
  !!  descr       - Print a more detailed description
  !!
  !!  cseti, csetr, csetc  - Set internal parameters, if any
  !
  type amg_s_base_aggregator_type
    ! Do we want to purge explicit zeros when aggregating?
    logical :: do_clean_zeros
  contains
    procedure, pass(ag) :: bld_tprol   => amg_s_base_aggregator_build_tprol
    procedure, pass(ag) :: mat_bld     => amg_s_base_aggregator_mat_bld
    procedure, pass(ag) :: mat_asb     => amg_s_base_aggregator_mat_asb
    procedure, pass(ag) :: bld_map     => amg_s_base_aggregator_bld_map
    procedure, pass(ag) :: update_next => amg_s_base_aggregator_update_next
    procedure, pass(ag) :: clone       => amg_s_base_aggregator_clone
    procedure, pass(ag) :: free        => amg_s_base_aggregator_free
    procedure, pass(ag) :: default     => amg_s_base_aggregator_default
    procedure, pass(ag) :: descr       => amg_s_base_aggregator_descr
    procedure, pass(ag) :: sizeof      => amg_s_base_aggregator_sizeof
    procedure, pass(ag) :: set_aggr_type => amg_s_base_aggregator_set_aggr_type
    procedure, nopass   :: fmt         => amg_s_base_aggregator_fmt
    procedure, pass(ag) :: cseti       => amg_s_base_aggregator_cseti
    procedure, pass(ag) :: csetr       => amg_s_base_aggregator_csetr
    procedure, pass(ag) :: csetc       => amg_s_base_aggregator_csetc
    generic, public     :: set         => cseti, csetr, csetc
    procedure, nopass   :: xt_desc     => amg_s_base_aggregator_xt_desc
  end type amg_s_base_aggregator_type

  abstract interface
    subroutine amg_s_soc_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
      import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_
      implicit none
      integer(psb_ipk_), intent(in)     :: iorder
      logical, intent(in)               :: clean_zeros
      type(psb_sspmat_type), intent(in) :: a
      type(psb_desc_type), intent(in)    :: desc_a
      real(psb_spk_), intent(in)         :: theta
      integer(psb_lpk_), allocatable, intent(out)  :: ilaggr(:),nlaggr(:)
      integer(psb_ipk_), intent(out)               :: info
    end subroutine amg_s_soc_map_bld
  end interface

  interface amg_ptap_bld
    subroutine amg_s_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
         & coo_prol,desc_cprol,coo_restr,info,desc_ax)
      import :: psb_s_csr_sparse_mat, psb_sspmat_type, psb_desc_type, &
           & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
      implicit none
      type(psb_s_csr_sparse_mat), intent(inout) :: a_csr
      type(psb_desc_type), intent(inout)         :: desc_a
      integer(psb_lpk_), intent(inout)           :: nlaggr(:)
      type(amg_sml_parms), intent(inout)         :: parms
      type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
      type(psb_desc_type), intent(inout)         :: desc_cprol
      type(psb_sspmat_type), intent(out)        :: ac
      integer(psb_ipk_), intent(out)             :: info
      type(psb_desc_type), intent(inout), optional :: desc_ax
    end subroutine amg_s_ptap_bld
  end interface amg_ptap_bld

  interface amg_rap
    subroutine amg_s_rap(a_csr,desc_a,nlaggr,parms,ac,&
         & coo_prol,desc_cprol,coo_restr,info)
      import :: psb_s_csr_sparse_mat, psb_sspmat_type, psb_desc_type, &
           & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
      implicit none
      type(psb_s_csr_sparse_mat), intent(inout) :: a_csr
      type(psb_desc_type), intent(inout)         :: desc_a
      integer(psb_lpk_), intent(inout)           :: nlaggr(:)
      type(amg_sml_parms), intent(inout)         :: parms
      type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
      type(psb_desc_type), intent(inout)         :: desc_cprol
      type(psb_sspmat_type), intent(out)        :: ac
      integer(psb_ipk_), intent(out)             :: info
    end subroutine amg_s_rap
  end interface amg_rap

contains

  subroutine amg_s_base_aggregator_cseti(ag,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_s_base_aggregator_type), intent(inout) :: ag
    character(len=*), intent(in)                  :: what
    integer(psb_ipk_), intent(in)                 :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    ! Do nothing
    info = 0
  end subroutine amg_s_base_aggregator_cseti

  subroutine amg_s_base_aggregator_csetr(ag,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_s_base_aggregator_type), intent(inout) :: ag
    character(len=*), intent(in)                  :: what
    real(psb_spk_), intent(in)                 :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    ! Do nothing
    info = 0
  end subroutine amg_s_base_aggregator_csetr

  subroutine amg_s_base_aggregator_csetc(ag,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_s_base_aggregator_type), intent(inout) :: ag
    character(len=*), intent(in)                  :: what
    character(len=*), intent(in)                 :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    ! Set clean zeros, or do nothing.
    select case (psb_toupper(trim(what)))
    case('AGGR_CLEAN_ZEROS')
      select case (psb_toupper(trim(val)))
      case('TRUE','T')
        ag%do_clean_zeros = .true.
      case('FALSE','F')
        ag%do_clean_zeros = .false.
      end select
    end select
    info = 0
  end subroutine amg_s_base_aggregator_csetc

  subroutine  amg_s_base_aggregator_update_next(ag,agnext,info)
    implicit none
    class(amg_s_base_aggregator_type), target, intent(inout) :: ag, agnext
    integer(psb_ipk_), intent(out)       :: info

    !
    ! Base version does nothing.
    !
    info = 0
  end subroutine amg_s_base_aggregator_update_next

  subroutine  amg_s_base_aggregator_clone(ag,agnext,info)
    implicit none
    class(amg_s_base_aggregator_type), intent(inout) :: ag
    class(amg_s_base_aggregator_type), allocatable, intent(inout) :: agnext
    integer(psb_ipk_), intent(out)       :: info

    info = 0
    if (allocated(agnext)) then
      call agnext%free(info)
      if (info == 0) deallocate(agnext,stat=info)
    end if
    if (info /= 0) return
    allocate(agnext,source=ag,stat=info)

  end subroutine amg_s_base_aggregator_clone

  subroutine  amg_s_base_aggregator_free(ag,info)
    implicit none
    class(amg_s_base_aggregator_type), intent(inout) :: ag
    integer(psb_ipk_), intent(out)       :: info

    info = psb_success_
    return
  end subroutine amg_s_base_aggregator_free

  subroutine  amg_s_base_aggregator_default(ag)
    implicit none
    class(amg_s_base_aggregator_type), intent(inout) :: ag
    ! Only one default setting
    ag%do_clean_zeros = .true.

    return
  end subroutine amg_s_base_aggregator_default

  function amg_s_base_aggregator_fmt() result(val)
    implicit none
    character(len=32)  :: val

    val = "Default aggregator "
  end function amg_s_base_aggregator_fmt

  function amg_s_base_aggregator_sizeof(ag) result(val)
    implicit none
    class(amg_s_base_aggregator_type), intent(in)  :: ag
    integer(psb_epk_)  :: val

    val = 1
  end function amg_s_base_aggregator_sizeof

  function amg_s_base_aggregator_xt_desc() result(val)
    implicit none
    logical  :: val

    val = .false.
  end function amg_s_base_aggregator_xt_desc

  subroutine  amg_s_base_aggregator_descr(ag,parms,iout,info)
    implicit none
    class(amg_s_base_aggregator_type), intent(in) :: ag
    type(amg_sml_parms), intent(in)   :: parms
    integer(psb_ipk_), intent(in)  :: iout
    integer(psb_ipk_), intent(out) :: info

    write(iout,*) 'Aggregator object type: ',ag%fmt()
    call parms%mldescr(iout,info)

    return
  end subroutine amg_s_base_aggregator_descr

  subroutine  amg_s_base_aggregator_set_aggr_type(ag,parms,info)
    implicit none
    class(amg_s_base_aggregator_type), intent(inout) :: ag
    type(amg_sml_parms), intent(in)   :: parms
    integer(psb_ipk_), intent(out) :: info

    ! Do nothing

    return
  end subroutine amg_s_base_aggregator_set_aggr_type

  !
  !> Function  bld_tprol:
  !! \memberof  amg_s_base_aggregator_type
  !! \brief  Build a tentative prolongator.
  !!         The routine will map the local matrix entries to aggregates.
  !!         The mapping is store in ILAGGR; for each local row index I,
  !!         ILAGGR(I) contains the index of the aggregate to which index I
  !!         will contribute, in global numbering.
  !!         Many aggregations produce a binary tentative prolongator, but some
  !!         do not, hence we also need the OP_PROL output.
  !!         AG_DATA is passed here just in case some of the
  !!         aggregators need it internally, most of them will ignore.
  !!
  !!  \param ag      The input aggregator object
  !!  \param parms   The auxiliary parameters object
  !!  \param ag_data Auxiliary global aggregation info
  !!  \param a       The local matrix part
  !!  \param desc_a  The descriptor
  !!  \param ilaggr  Output aggregation map
  !!  \param nlaggr  Sizes of ilaggr on all processes
  !!  \param op_prol The tentative prolongator operator
  !!  \param info    Return code
  !!
  !
  subroutine  amg_s_base_aggregator_build_tprol(ag,parms,ag_data,&
       & a,desc_a,ilaggr,nlaggr,t_prol,info)
    use psb_base_mod
    implicit none
    class(amg_s_base_aggregator_type), target, intent(inout) :: ag
    type(amg_sml_parms), intent(inout)  :: parms
    type(amg_saggr_data), intent(in)    :: ag_data
    type(psb_sspmat_type), intent(inout) :: a
    type(psb_desc_type), intent(inout)     :: desc_a
    integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
    type(psb_lsspmat_type), intent(out)  :: t_prol
    integer(psb_ipk_), intent(out)      :: info

    integer(psb_ipk_) :: err_act
    character(len=20) :: name='s_base_aggregator_build_tprol'

    call psb_erractionsave(err_act)

    info = psb_err_missing_override_method_
    call psb_errpush(info,name)
    goto 9999

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return

  end subroutine amg_s_base_aggregator_build_tprol

  !
  !> Function   mat_bld
  !! \memberof  amg_s_base_aggregator_type
  !! \brief     Build prolongator/restrictor/coarse matrix.
  !!
  !!
  !!  \param ag        The input aggregator object
  !!  \param parms     The auxiliary parameters object
  !!  \param a         The local matrix part
  !!  \param desc_a    The descriptor
  !!  \param ilaggr    Aggregation map
  !!  \param nlaggr    Sizes of ilaggr on all processes
  !!  \param ac        On output the coarse matrix
  !!  \param op_prol   On input, the  tentative prolongator operator, on output
  !!                   the final prolongator
  !!  \param op_restr  On output, the restrictor operator;
  !!                   in many cases it is the transpose of the prolongator.
  !!  \param info    Return code
  !!
  subroutine  amg_s_base_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
       & ac,desc_ac,op_prol,op_restr,t_prol,info)
    use psb_base_mod
    implicit none
    class(amg_s_base_aggregator_type), target, intent(inout) :: ag
    type(amg_sml_parms), intent(inout)   :: parms
    type(psb_sspmat_type), intent(in)    :: a
    type(psb_desc_type), intent(inout)   :: desc_a
    integer(psb_lpk_), intent(inout)     :: ilaggr(:), nlaggr(:)
    type(psb_lsspmat_type), intent(inout) :: t_prol
    type(psb_sspmat_type), intent(out)   :: op_prol, ac,op_restr
    type(psb_desc_type), intent(inout)   :: desc_ac
    integer(psb_ipk_), intent(out)       :: info
    integer(psb_ipk_) :: err_act
    character(len=20) :: name='s_base_aggregator_mat_bld'

    call psb_erractionsave(err_act)

    info = psb_err_missing_override_method_
    call psb_errpush(info,name)
    goto 9999

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return
  end subroutine amg_s_base_aggregator_mat_bld

  !
  !> Function   mat_asb
  !! \memberof  amg_s_base_aggregator_type
  !! \brief     Build prolongator/restrictor/coarse matrix.
  !!
  !!
  !!  \param ag        The input aggregator object
  !!  \param parms     The auxiliary parameters object
  !!  \param a         The local matrix part
  !!  \param desc_a    The descriptor
  !!  \param ilaggr    Aggregation map
  !!  \param nlaggr    Sizes of ilaggr on all processes
  !!  \param ac        On output the coarse matrix
  !!  \param op_prol   On input, the  tentative prolongator operator, on output
  !!                   the final prolongator
  !!  \param op_restr  On output, the restrictor operator;
  !!                   in many cases it is the transpose of the prolongator.
  !!  \param info    Return code
  !!
  subroutine  amg_s_base_aggregator_mat_asb(ag,parms,a,desc_a,&
       & ac,desc_ac, op_prol,op_restr,info)
    use psb_base_mod
    implicit none
    class(amg_s_base_aggregator_type), target, intent(inout) :: ag
    type(amg_sml_parms), intent(inout)   :: parms
    type(psb_sspmat_type), intent(in)    :: a
    type(psb_desc_type), intent(inout)     :: desc_a
    type(psb_sspmat_type), intent(inout) :: op_prol,ac,op_restr
    type(psb_desc_type), intent(inout)      :: desc_ac
    integer(psb_ipk_), intent(out)       :: info
    integer(psb_ipk_) :: err_act
    character(len=20) :: name='s_base_aggregator_mat_asb'

    call psb_erractionsave(err_act)

    info = psb_err_missing_override_method_
    call psb_errpush(info,name)
    goto 9999

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return
  end subroutine amg_s_base_aggregator_mat_asb

  !
  !> Function   bld_map
  !! \memberof  amg_s_base_aggregator_type
  !! \brief     Build linear map between hierarchy levels
  !!
  !!
  !!  \param ag        The input aggregator object
  !!  \param desc_a    The fine space descriptor
  !!  \param desc_ac   The coarse space descriptor
  !!  \param ilaggr    Aggregation map vector
  !!  \param nlaggr    Sizes of ilaggr on all processes
  !!  \param op_prol   The  prolongator operator
  !!  \param op_restr  The restrictor operator
  !!  \param map       The output map
  !!  \param info    Return code
  !!
  subroutine  amg_s_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
       & op_restr,op_prol,map,info)
    use psb_base_mod
    implicit none
    class(amg_s_base_aggregator_type), target, intent(inout) :: ag
    type(psb_desc_type), intent(in), target :: desc_a, desc_ac
    integer(psb_lpk_), intent(inout)     :: ilaggr(:), nlaggr(:)
    type(psb_sspmat_type), intent(inout)   :: op_restr, op_prol
    type(psb_slinmap_type), intent(out)    :: map
    integer(psb_ipk_), intent(out)       :: info
    integer(psb_ipk_) :: err_act
    character(len=20) :: name='s_base_aggregator_bld_map'

    call psb_erractionsave(err_act)
    !
    ! Copy the prolongation/restriction matrices into the descriptor map.
    !  op_restr => PR^T   i.e. restriction  operator
    !  op_prol => PR     i.e. prolongation operator
    !
    !  WARNING: need to check whether the  copy into IOP_RESTR/IOP_PROL
    !  is safe or not.
    !
    !  This default implementation reuses desc_a/desc_ac through
    !  pointers in the map structure.
    !
    map = psb_linmap(psb_map_aggr_,desc_a,&
         & desc_ac,op_restr,op_prol,ilaggr,nlaggr)

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return
  end subroutine amg_s_base_aggregator_bld_map

end module amg_s_base_aggregator_mod


text:
cc      MODULE UELUTIL
cc      CONTAINS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C11111111122222222223333333333444444444455555555556666666666777777777777
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C       --U S E R   E L E M E N T    S U B R O U T I N E S---          C
C                                                                      C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C  UELEXP9.for                                                         C
C  9-NODED ISOPARAMETRIC ELEMENT WITH FULL GAUSS INTEGRATION           C
C  USER ELEMENT SUBROUTINE FOR EXPLICIT DYNAMIC ANALYSIS               C
C                                                                      C
C  GIVEN THE DISPLACEMENT AT TIME T AND THE DISPLACEMENT INCREMENT     C
C  IT RETURNS IN THE RHS() VECTOR THE INERTIAL FORCE CONTRIBUTION      C
C  AT TIME T+DT AND IN THE ANMASS() VECTOR THE DIAGONALIZED MASS       C
C  MATRIX                                                              C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UELEXP9(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,
     1                   AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,V,A,JELEM,
     2                   DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,NGPTS=9,
     1           THREE=3.D0, NDI=3)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),A(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     3          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL),UP(NDOFEL),
     4          DDSDDE(NDI,NDI),B(NDI,NDOFEL),BT(NDOFEL,NDI),
     5          FRST2(NDOFEL,NDOFEL),FRST1(NDOFEL,NDOFEL),XP(2,NGPTS),
     6          XW(NGPTS),AUX1(NDI,NDOFEL)

C     Mass density.

      RO=PROPS(3)

C     Rayleigh damping parameters.

      AALFA=PROPS(1)
      ABETA=PROPS(2)

      RHS=0.0D0
      ANMASS=0.0D0
	  DDSDDE=0.0D0
	  B =0.0D0
	  BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT(DDSDDE,NDI,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS3X3(XP,XW)
      NGPT=9
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        ALF=XW(NN)

C       Computes B matrix.

        CALL STDM9(NNODE,NDOFEL,NDI,COORDS,MCRD,B,DDET,RII,SII,
     1             XBAR)

        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

	  	AUX1=MATMUL(DDSDDE,B)
	  	BT=TRANSPOSE(B)
	  	FRST1=MATMUL(BT,AUX1)
C
        CALL AMASS9(FRST2,NDOFEL,RII,SII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrix.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes elastic,inertial and damping nodal forces.

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

      RETURN

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C  UELEXP8.for                                                         C
C  8-NODED ISOPARAMETRIC ELEMENT WITH FULL GAUSS INTEGRATION           C
C  USER ELEMENT SUBROUTINE FOR EXPLICIT DYNAMIC ANALYSIS               C
C                                                                      C
C  GIVEN THE DISPLACEMENT AT TIME T AND THE DISPLACEMENT INCREMENT     C
C  IT RETURNS IN THE RHS() VECTOR THE INERTIAL FORCE CONTRIBUTION      C
C  AT TIME T+DT AND IN THE ANMASS() VECTOR THE DIAGONALIZED MASS       C
C  MATRIX                                                              C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UELEXP8(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,
     1                   AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,V,A,JELEM,
     2                   DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,NGPTS=9,
     1           THREE=3.D0, NDI=3)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),A(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     3          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL),UP(NDOFEL),
     4          DDSDDE(NDI,NDI),B(NDI,NDOFEL),BT(NDOFEL,NDI),
     5          FRST2(NDOFEL,NDOFEL),FRST1(NDOFEL,NDOFEL),XP(2,NGPTS),
     6          XW(NGPTS),AUX1(NDI,NDOFEL)

C     Mass density.

      RO=PROPS(3)

C     Rayleigh damping parameters.

      AALFA=PROPS(1)
      ABETA=PROPS(2)

      RHS=0.0D0
      ANMASS=0.0D0
	  DDSDDE=0.0D0
	  B =0.0D0
	  BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT(DDSDDE,NDI,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS3X3(XP,XW)
      NGPT=9
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        ALF=XW(NN)

C       Computes B matrix.

        CALL STDM8(NNODE,NDOFEL,NDI,COORDS,MCRD,B,DDET,RII,SII,
     1             XBAR)
        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

	  	AUX1=MATMUL(DDSDDE,B)
	  	BT=TRANSPOSE(B)
	  	FRST1=MATMUL(BT,AUX1)
C
        CALL AMASS8(FRST2,NDOFEL,RII,SII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrix.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes elastic,inertial and damping nodal forces.

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

      RETURN

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C  UELEXP6.for                                                         C
C  6-NODED ISOPARAMETRIC ELEMENT WITH FULL GAUSS INTEGRATION           C
C  USER ELEMENT SUBROUTINE FOR EXPLICIT DYNAMIC ANALYSIS               C
C                                                                      C
C  GIVEN THE DISPLACEMENT AT TIME T AND THE DISPLACEMENT INCREMENT     C
C  IT RETURNS IN THE RHS() VECTOR THE INERTIAL FORCE CONTRIBUTION      C
C  AT TIME T+DT AND IN THE ANMASS() VECTOR THE DIAGONALIZED MASS       C
C  MATRIX                                                              C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UELEXP6(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,
     1                   AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,V,A,JELEM,
     2                   DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,NGPTS=7,
     1           THREE=3.D0, NDI=3)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),A(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     3          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL),UP(NDOFEL),
     4          DDSDDE(NDI,NDI),B(NDI,NDOFEL),BT(NDOFEL,NDI),
     5          FRST2(NDOFEL,NDOFEL),FRST1(NDOFEL,NDOFEL),XP(2,NGPTS),
     6          XW(NGPTS),AUX1(NDI,NDOFEL)

C     Mass density.

      RO=PROPS(3)

C     Rayleigh damping parameters.

      AALFA=PROPS(1)
      ABETA=PROPS(2)

      RHS=0.0D0
      ANMASS=0.0D0
	  DDSDDE=0.0D0
	  B =0.0D0
	  BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT(DDSDDE,NDI,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS7(XP,XW)
      NGPT=7
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        ALF=XW(NN)

C       Computes B matrix.

        CALL STDM6(NNODE,NDOFEL,NDI,COORDS,MCRD,B,DDET,RII,SII,
     1             XBAR)
        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

	  	AUX1=MATMUL(DDSDDE,B)
	  	BT=TRANSPOSE(B)
	  	FRST1=MATMUL(BT,AUX1)
C
        CALL AMASS6(FRST2,NDOFEL,RII,SII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrix.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes elastic,inertial and damping nodal forces.

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C  UELEXP4.for                                                         C
C  4-NODED ISOPARAMETRIC ELEMENT WITH FULL GAUSS INTEGRATION           C
C  USER ELEMENT SUBROUTINE FOR EXPLICIT DYNAMIC ANALYSIS               C
C                                                                      C
C  GIVEN THE DISPLACEMENT AT TIME T AND THE DISPLACEMENT INCREMENT     C
C  IT RETURNS IN THE RHS() VECTOR THE INERTIAL FORCE CONTRIBUTION      C
C  AT TIME T+DT AND IN THE ANMASS() VECTOR THE DIAGONALIZED MASS       C
C  MATRIX                                                              C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UELEXP4(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,
     1                   AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,V,A,JELEM,
     2                   DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,NGPTS=4,
     1           THREE=3.D0, NDI=3)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),A(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     3          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL),UP(NDOFEL),
     4          DDSDDE(NDI,NDI),B(NDI,NDOFEL),BT(NDOFEL,NDI),
     5          FRST2(NDOFEL,NDOFEL),FRST1(NDOFEL,NDOFEL),XP(2,NGPTS),
     6          XW(NGPTS),AUX1(NDI,NDOFEL)

C     Mass density.

      RO=PROPS(3)
C
C     Rayleigh damping parameters.

      AALFA=PROPS(1)
      ABETA=PROPS(2)

      RHS=0.0D0
      ANMASS=0.0D0
	  DDSDDE=0.0D0
	  B =0.0D0
	  BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT(DDSDDE,NDI,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS2X2(XP,XW)
      NGPT=4
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        ALF=XW(NN)

C       Computes B matrix.

        CALL STDM4(NNODE,NDOFEL,NDI,COORDS,MCRD,B,DDET,RII,SII,
     1             XBAR)
        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

	  	AUX1=MATMUL(DDSDDE,B)
	  	BT=TRANSPOSE(B)
	  	FRST1=MATMUL(BT,AUX1)

        CALL AMASS4(FRST2,NDOFEL,RII,SII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrix.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes elastic,inertial and damping nodal forces.

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE UEL9INCOT(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,  C
C    1                   COORDS,MCRD,NNODE,U,DU,V,JELEM,DT)            C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UEL9INCOT(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                  NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,V,A,
     2                  JELEM,DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,TWO=2.D0,NGPTS=9,
     1           THREE=3.D0,NDI=3,ONE=1.0D0)

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),AWAVE(5),
     1          COORDS(MCRD,NNODE),U(NDOFEL),DU(NDOFEL),V(NDOFEL),
     2          UP(NDOFEL),DDSDDE(NDI,NDI),B(NDI,NDOFEL),
     3          BT(NDOFEL,NDI),FRST2(NDOFEL,NDOFEL),A(NDOFEL),
     4          FRST1(NDOFEL,NDOFEL),XP(2,NGPTS),XW(NGPTS),
     5          AUX1(NDI,NDOFEL),R00(NDOFEL),RELA(NDOFEL),
     6          RINE(NDOFEL),RDAM(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     7          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL)

      DIMENSION IPROPS(NIPROPS)

C     Mass density

      RO=PROPS(3)

C     Rayleigh damping constants.

      AALFA=PROPS(1)
      ABETA=PROPS(2)
      IDF=IPROPS(1)

C     Clears arrays.

      RHS=0.0D0
      ANMASS=0.0D0
      DDSDDE=0.0D0
      B=0.0D0
      BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      R00=0.0D0
      RELA=0.0D0
      RINE=0.0D0
      RDAM=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT(DDSDDE,NDI,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS3X3(XP,XW)
      NGPT=9
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        ALF=XW(NN  )

C       Computes B and mass matrix.

        CALL STDM9(NNODE,NDOFEL,NDI,COORDS,MCRD,B,DDET,RII,SII,
     1             XBAR)
        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

	  	AUX1=MATMUL(DDSDDE,B)
	  	BT=TRANSPOSE(B)
	  	FRST1=MATMUL(BT,AUX1)

        CALL AMASS9(FRST2,NDOFEL,RII,SII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrices.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes the contribution -Kij*Uj-MII*Ai-alpha*Kij*Vj.....

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

C     Computes the incoming wave field and mekes it into effective
C     forces R00(n).

      CALL INCOMINGDIS9(IDF,AMATRX,COORDS,PROPS,NPROPS,MCRD,AWAVE,
     1                 IWAVE,NNODE,NDOFEL,R00,DT,KINC,NINCR,IOUT)

C     Updates the effective loads
C     loads vector RHS=R00(n)-Kij*Uj-MII*Ai-alpha*Kij*Vj....

      RHS=RHS+R00

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE UEL8INCOT(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,C
C    1                   COORDS,MCRD,NNODE,U,DU,V,JELEM,DT)            C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UEL8INCOT(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                  NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,V,A,
     2                  JELEM,DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,TWO=2.D0,NGPTS=9,
     1           THREE=3.D0,NDI=3,ONE=1.0D0)

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),AWAVE(5),
     1          COORDS(MCRD,NNODE),U(NDOFEL),DU(NDOFEL),V(NDOFEL),
     2          UP(NDOFEL),DDSDDE(NDI,NDI),B(NDI,NDOFEL),
     3          BT(NDOFEL,NDI),FRST2(NDOFEL,NDOFEL),A(NDOFEL),
     4          FRST1(NDOFEL,NDOFEL),XP(2,NGPTS),XW(NGPTS),
     5          AUX1(NDI,NDOFEL),R00(NDOFEL),RELA(NDOFEL),
     6          RINE(NDOFEL),RDAM(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     7          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL)

      DIMENSION IPROPS(NIPROPS)

C     Mass density

      RO=PROPS(3)

C     Rayleigh damping constants.

      AALFA=PROPS(1)
      ABETA=PROPS(2)
      IDF=IPROPS(1)

C     Clears arrays.

      RHS=0.0D0
      ANMASS=0.0D0
      DDSDDE=0.0D0
      B=0.0D0
      BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      R00=0.0D0
      RELA=0.0D0
      RINE=0.0D0
      RDAM=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT(DDSDDE,NDI,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS3X3(XP,XW)
      NGPT=9
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        ALF=XW(NN  )

C       Computes B and mass matrix.

        CALL STDM8(NNODE,NDOFEL,NDI,COORDS,MCRD,B,DDET,RII,SII,
     1             XBAR)
        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

	  	AUX1=MATMUL(DDSDDE,B)
	  	BT=TRANSPOSE(B)
	  	FRST1=MATMUL(BT,AUX1)

        CALL AMASS8(FRST2,NDOFEL,RII,SII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrices.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes the contribution -Kij*Uj-MII*Ai-alpha*Kij*Vj.....

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

C     Computes the incoming wave field and mekes it into effective
C     forces R00(n).

      CALL INCOMINGDIS8(IDF,AMATRX,COORDS,PROPS,NPROPS,MCRD,AWAVE,
     1                 IWAVE,NNODE,NDOFEL,R00,DT,KINC,NINCR,IOUT)

C     Updates the effective loads
C     loads vector RHS=R00(n)-Kij*Uj-MII*Ai-alpha*Kij*Vj....

      RHS=RHS+R00

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C   SUBROUTINE UEL4INCOT(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,C
C    1                   COORDS,MCRD,NNODE,U,DU,V,JELEM,DT)            C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UEL4INCOT(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                  NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,V,A,
     2                  JELEM,DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,NGPTS=4,
     1           THREE=3.D0, NDI=3)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),AWAVE(5),
     1          COORDS(MCRD,NNODE),U(NDOFEL),DU(NDOFEL),V(NDOFEL),
     2          UP(NDOFEL),DDSDDE(NDI,NDI),B(NDI,NDOFEL),
     3          BT(NDOFEL,NDI),FRST2(NDOFEL,NDOFEL),A(NDOFEL),
     4          FRST1(NDOFEL,NDOFEL),XP(2,NGPTS),XW(NGPTS),
     5          AUX1(NDI,NDOFEL),R00(NDOFEL),RELA(NDOFEL),
     6          RINE(NDOFEL),RDAM(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     7          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL)

      DIMENSION IPROPS(NIPROPS)

C     Mass density.

      RO=PROPS(3)
C
C     Rayleigh damping parameters.

      AALFA=PROPS(1)
      ABETA=PROPS(2)
      IDF=IPROPS(1)

C     Clears arrays.

      RHS=0.0D0
      ANMASS=0.0D0
      DDSDDE=0.0D0
      B=0.0D0
      BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      R00=0.0D0
      RELA=0.0D0
      RINE=0.0D0
      RDAM=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT(DDSDDE,NDI,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS2X2(XP,XW)
      NGPT=4
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        ALF=XW(NN)

C       Computes B matrix.

        CALL STDM4(JELEM,NNODE,NDOFEL,NDI,COORDS,MCRD,B,DDET,RII,SII,
     1             XBAR)
        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

	  	AUX1=MATMUL(DDSDDE,B)
	  	BT=TRANSPOSE(B)
	  	FRST1=MATMUL(BT,AUX1)

        CALL AMASS4(FRST2,NDOFEL,RII,SII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrix.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes elastic,inertial and damping nodal forces.

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

C     Computes the incoming wave field and mekes it into effective
C     forces R00(n).

      CALL INCOMINGDIS4(IDF,AMATRX,COORDS,PROPS,NPROPS,MCRD,AWAVE,
     1                  IWAVE,NNODE,NDOFEL,R00,DT,KINC,NINCR,IOUT)

C     Updates the effective loads
C     loads vector RHS=R00(n)-Kij*Uj-MII*Ai-alpha*Kij*Vj....

      RHS=RHS+R00

      RETURN

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C    SUBROUTINE UELEDASH3(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,C
C    1                    COORDS,MCRD,NNODE,U,DU,V,JELEM,DT)           C
C                                                                      C
C    ABSORBING DASHPOT ELEMENT                                         C
C                                                                      C
C                                                                      C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE UELEDASH3(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                    NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,
     2                    U,DU,V,A,JELEM,DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),AMATRX(NDOFEL,NDOFEL),DDSDDE(2,2),
     3          XGP(6),WGP(6),DDSDDET(2,2), DAUX(2,2),AUX(2,6),
     4          AUXT(6,2),RL(6,6),RLL(6,2),RM(2,2),RMT(2,2),A(NDOFEL)

C     CALL CLEARV(RHS,NDOFEL)
C     CALL CLEARV(ANMASS,NDOFEL)
C     CALL CLEAR(AMATRX,NDOFEL,NDOFEL)
      AMATRX=0.0D0

C     Wave propagation velocities

      AALFA=PROPS(1)
      ABETA=PROPS(2)
      RO=PROPS(3)

C     Constitutive tensor

      DDSDDE(1,2)=0.0D0
      DDSDDE(2,1)=0.0D0
      DDSDDE(1,1)=RO*AALFA
      DDSDDE(2,2)=RO*ABETA

      XI=COORDS(1,1)
      XM=COORDS(1,2)
      XJ=COORDS(1,3)
      YI=COORDS(2,1)
      YM=COORDS(2,2)
      YJ=COORDS(2,3)
      DXS=(XI-XJ)**2
      DYS=(YI-YJ)**2
      DL=DSQRT(DXS+DYS)
      DETJAC=DL/2.0

      CALL NORMVCT(XI,XJ,YI,YJ,CT,ST)
      RM(1,1)=  CT
      RM(1,2)=  ST
      RM(2,1)= -ST
      RM(2,2)=  CT

      RMT=TRANSPOSE(RM)
      DAUX=MATMUL(RMT,DDSDDE)
      DDSDDET=MATMUL(DAUX,RM)

      CALL GAUSS1D(XGP,WGP)

      DO N=1,6
        ETA =XGP(N)
        ALF =WGP(N)
        CALL SHAPE1D(AUXT,ETA)

        RLL=MATMUL(AUXT, DDSDDET)
	  	AUX=TRANSPOSE(AUXT)
	  	RL=MATMUL(RLL,AUX)
        AMATRX=AMATRX-RL*ALF*DETJAC
      END DO

      RHS=MATMUL(AMATRX,V)

      RETURN

      END

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C    SUBROUTINE UELDASH3D(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPSC
C    1                    ,COORDS,MCRD,NNODE,U,DU,V,JELEM,DT)          C
C                                                                      C
C    ABSORBING 3D DASHPOT ELEMENT                                      C
C                                                                      C
C                                                                      C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE UELDASH3D(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                    NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,
     2                    U,DU,V,A,JELEM,DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,
     1           NGPTS=4)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),AMATRX(NDOFEL,NDOFEL),DDSDDE(3,3),
     3          XGP(2,NGPTS),WGP(NGPTS),DDSDDET(3,3),DAUX(2,2),
     4          AUX(3,12),AUXT(12,3),RL(12,12),RLL(12,3),RM(3,3),
     5			RMT(3,3),A(NDOFEL),P1(3), P2(3), P3(3), P4(3),
     6			V1(3),V2(3),V3(3),V4(3),P1T(3),P2T(3),P3T(3),P4T(3),
     7			XY(2,4), P(2,4), XJ(2,2)

      AMATRX=0.0D0

C     Wave propagation velocities

      AALFA=PROPS(1)
      ABETA=PROPS(2)
      RO=PROPS(3)

C     Constitutive tensor

      DDSDDE=0.0D0

      DDSDDE(3,3)=RO*ABETA
      DDSDDE(2,2)=RO*ABETA
	  DDSDDE(1,1)=RO*AALFA

	  CALL TRANSFORMA(COORDS, RM, XY)

      RMT=TRANSPOSE(RM)

      DDSDDET= MATMUL(RM, DDSDDE)

	  CALL GPOINTS2X2(XGP,WGP)
	  DDSDDE = MATMUL(DDSDDET, RMT)

	  SMT= 0.0D0

      DO N=1,NGPTS

      	RII= XGP(1,N)
      	SII= XGP(2,N)
        ALF= WGP(N)

        CALL SFDER4(1, NDOFEL, 4, RII, SII, P)
        CALL JACOPER(4,XJ,XX,P)
	    DETJAC= (XJ(1,1)*XJ(2,2)-XJ(1,2)*XJ(2,1))
	    CALL SHAPEF4DASH(AUXT, RII, SII, KINC)

	  	RLL= MATMUL(AUXT, DDSDDE)
	  	AUX= TRANSPOSE(AUXT)
	  	RL= MATMUL(RLL, AUX)
        AMATRX=AMATRX-RL*ALF*DETJAC
      END DO

      CALL MAVEC(AMATRX,NDOFEL,NDOFEL,V,RHS)

      RETURN

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C    SUBROUTINE UELEDASHINC(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,C
C    1                    COORDS,MCRD,NNODE,U,DU,V,JELEM,DT)           C
C                                                                      C
C    TRANSMITING DASHPOT ELEMENT USED TO EXCITATE THE BASE OF          C
C    A MODEL WITH A SV WAVE WITH VERTICAL INCIDENCE                    C
C                                                                      C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE UELDASHINC(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                      NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,
     2                      U,DU,V,A,JELEM,DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),
     2          U(NDOFEL),DU(NDOFEL),V(NDOFEL),
     3          AMATRX(NDOFEL,NDOFEL),DDSDDE(2,2),XGP(6),VIN(6),
     4          WGP(6),AUX(2,6),AUXT(6,2),RL(6,6),RLL(6,2),A(NDOFEL)

      AMATRX=0.0D0
      VIN=0.0D0
      PPV=0.0D0

C     Wave propagation velocities

      AALFA=PROPS(1)
      ABETA=PROPS(2)
      RO=PROPS(3)

C     Constitutive tensor

      DDSDDE(1,2)=0.0D0
      DDSDDE(2,1)=0.0D0
      DDSDDE(1,1)=RO*ABETA
      DDSDDE(2,2)=0.0D0

      XI=COORDS(1,1)
      XM=COORDS(1,2)
      XJ=COORDS(1,3)
      YI=COORDS(2,1)
      YM=COORDS(2,2)
      YJ=COORDS(2,3)
      DXS=(XI-XJ)**2
      DYS=(YI-YJ)**2
      DL=DSQRT(DXS+DYS)
      DETJAC=DL/2.0D0

      CALL GAUSS1D(XGP,WGP)

      DO N=1,6
        ETA =XGP(N)
        ALF=WGP(N)
        CALL SHAPE1D(AUXT,ETA)

        RLL=MATMUL(AUXT,DDSDDE)
        AUX=TRANSPOSE(AUXT)
        RL=MATMUL(RLL,AUX)

        AMATRX=AMATRX-RL*ALF*DETJAC
      END DO

      CALL UPULSEV(PPV,NINCR,DT,AWAVE,KINC)

      VIN(1)=PPV
      VIN(3)=PPV
      VIN(5)=PPV

      RHS=MATMUL(AMATRX,VIN)

      RETURN

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C    SUBROUTINE UELEDASH2(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,NIPROPS,C
C    1                    COORDS,MCRD,NNODE,U,DU,V,JELEM,DT)           C
C                                                                      C
C    ABSORBING DASHPOT ELEMENT                                         C
C                                                                      C
C                                                                      C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE UELEDASH2(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                    NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,
     2                    U,DU,V,A,JELEM,DT,KINC,NINCR)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),AMATRX(NDOFEL,NDOFEL),DDSDDE(2,2),
     3          XGP(4),WGP(4),DDSDDET(2,2), DAUX(2,2),AUX(2,4),
     4          AUXT(4,2),RL(4,4),RLL(4,2),RM(2,2),RMT(2,2),A(NDOFEL)

      AMATRX=0.0D0

C     Wave propagation velocities

      AALFA=PROPS(1)
      ABETA=PROPS(2)
      RO=PROPS(3)

C     Constitutive tensor

      DDSDDE(1,2)=0.0D0
      DDSDDE(2,1)=0.0D0
      DDSDDE(1,1)=RO*AALFA
      DDSDDE(2,2)=RO*ABETA

      XI=COORDS(1,1)
      XM=COORDS(1,2)
      XJ=COORDS(1,3)
      YI=COORDS(2,1)
      YM=COORDS(2,2)
      YJ=COORDS(2,3)
      DXS=(XI-XJ)**2
      DYS=(YI-YJ)**2
      DL=DSQRT(DXS+DYS)
      DETJAC=DL/2.0

      CALL NORMVCT(XI,XJ,YI,YJ,CT,ST)
      RM(1,1)=  CT
      RM(1,2)=  ST
      RM(2,1)= -ST
      RM(2,2)=  CT

      RMT=TRANSPOSE(RM)
      DAUX=MATMUL(RMT,DDSDDE)
      DDSDDET=MATMUL(DAUX,RM)

      CALL GAUSS1D(XGP,WGP)

      DO N=1,4
        ETA =XGP(N)
        ALF=WGP(N)
        CALL SHAPE1D(AUXT,ETA)

        RLL=MATMUL(AUXT,DDSDDET)
        AUX=TRANSPOSE(AUXT)
        RL=MATMUL(RLL,AUX)

        AMATRX=AMATRX-RL*ALF*DETJAC
      END DO

      RHS=MATMUL(AMATRX,V)

      RETURN

      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C  UEL3DEXP8.for                                                       C
C  8-NODED ISOPARAMETRIC 3D ELEMENT WITH FULL GAUSS INTEGRATION        C
C  USER ELEMENT SUBROUTINE FOR EXPLICIT DYNAMIC ANALYSIS               C
C                                                                      C
C  GIVEN THE DISPLACEMENT AT TIME T AND THE DISPLACEMENT INCREMENT     C
C  IT RETURNS IN THE RHS() VECTOR THE INERTIAL FORCE CONTRIBUTION      C
C  AT TIME T+DT AND IN THE ANMASS() VECTOR THE DIAGONALIZED MASS       C
C  MATRIX                                                              C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UEL3DEXP8(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                   NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,
     2                   V,A,JELEM,DT,KINC,NINCR,INOD)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,NTENS=6,TWO=2.D0,
     1           NGPTS=8,THREE=3.D0, NDI=3)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),A(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     3          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL),UP(NDOFEL),
     4          DDSDDE(NTENS,NTENS),B(NTENS,NDOFEL),BT(NDOFEL,NTENS),
     5          FRST2(NDOFEL,NDOFEL),FRST1(NDOFEL,NDOFEL),XP(3,NGPTS),
     6          XW(NGPTS),AUX1(NTENS,NDOFEL)

C     Mass density.

      RO=PROPS(3)

C     Rayleigh damping parameters.

      AALFA=PROPS(1)
      ABETA=PROPS(2)

      RHS=0.0D0
      ANMASS=0.0D0
      DDSDDE=0.0D0
      B=0.0D0
      BT=0.0D0
      FRST2=0.0D0
      FRST1=0.0D0
      XP=0.0D0
      XW=0.0D0
      AUX1=0.0D0
      AMATRX=0.0D0
      AMASS=0.0D0
      CDMAT=0.0D0

C       Elasticity tensor.

	  CALL UMAT3D(DDSDDE,NDI,NTENS,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTS3D2X2(XP,XW)
      NGPT=8
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        TII=XP(3,NN)
        ALF=XW(NN)

C       Computes B matrix.

        CALL STDM3D8(NNODE,NDOFEL,NTENS,COORDS,MCRD,B,DDET,RII,SII,
     1               TII,XBAR)

        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

        AUX1=MATMUL(DDSDDE,B)
        BT=TRANSPOSE(B)
        FRST1=MATMUL(BT,AUX1)
C
        CALL AMASS3D8(FRST2,NDOFEL,RII,SII,TII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrix.

        FRST1=FRST1*(ALF*DDET*XBAR)
        FRST2=FRST2*(ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        AMATRX=AMATRX+FRST1
        AMASS=AMASS+FRST2

      END DO

C     Computes elastic,inertial and damping nodal forces.

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

      RETURN

      END
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C  UEL3DEXP8.for                                                       C
C  8-NODED ISOPARAMETRIC 3D ELEMENT WITH FULL GAUSS INTEGRATION        C
C  USER ELEMENT SUBROUTINE FOR EXPLICIT DYNAMIC ANALYSIS               C
C                                                                      C
C  GIVEN THE DISPLACEMENT AT TIME T AND THE DISPLACEMENT INCREMENT     C
C  IT RETURNS IN THE RHS() VECTOR THE INERTIAL FORCE CONTRIBUTION      C
C  AT TIME T+DT AND IN THE ANMASS() VECTOR THE DIAGONALIZED MASS       C
C  MATRIX                                                              C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE UEL3DWED10(RHS,ANMASS,NDOFEL,PROPS,NPROPS,IPROPS,
     1                   NIPROPS,AWAVE,IWAVE,COORDS,MCRD,NNODE,U,DU,
     2                   V,A,JELEM,DT,KINC,NINCR,INOD)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ZERO=0.D0,HALF=0.5D0,ONE=1.D0,NTENS=6,TWO=2.D0,
     1           NGPTS=8,THREE=3.D0, NDI=3)

C     Parameter arrays from UEL.f

      DIMENSION RHS(NDOFEL),ANMASS(NDOFEL),PROPS(NPROPS),
     1          IPROPS(NIPROPS),AWAVE(5),COORDS(MCRD,NNODE),U(NDOFEL),
     2          DU(NDOFEL),V(NDOFEL),A(NDOFEL),AMATRX(NDOFEL,NDOFEL),
     3          AMASS(NDOFEL,NDOFEL),CDMAT(NDOFEL,NDOFEL),UP(NDOFEL),
     4          DDSDDE(NTENS,NTENS),B(NTENS,NDOFEL),BT(NDOFEL,NTENS),
     5          FRST2(NDOFEL,NDOFEL),FRST1(NDOFEL,NDOFEL),XP(3,NGPTS),
     6          XW(NGPTS),AUX1(NTENS,NDOFEL)

C     Mass density.

      RO=PROPS(3)

C     Rayleigh damping parameters.

      AALFA=PROPS(1)
      ABETA=PROPS(2)

      CALL CLEARV(RHS,NDOFEL)
      CALL CLEARV(ANMASS,NDOFEL)
      CALL CLEAR(DDSDDE,NTENS,NTENS)
      CALL CLEAR(B,NTENS,NDOFEL)
      CALL CLEAR(BT,NDOFEL,NTENS)
      CALL CLEAR(FRST2,NDOFEL,NDOFEL)
      CALL CLEAR(FRST1,NDOFEL,NDOFEL)
      CALL CLEAR(XP,3,NGPTS)
      CALL CLEARV(XW,NGPTS)
      CALL CLEAR(AUX1,NTENS,NDOFEL)
      CALL CLEAR(AMATRX,NDOFEL,NDOFEL)
      CALL CLEAR(AMASS,NDOFEL,NDOFEL)
      CALL CLEAR(CDMAT,NDOFEL,NDOFEL)

C       Elasticity tensor.

	  CALL UMAT3D(DDSDDE,NDI,NTENS,PROPS,NPROPS)

C     Generates Gauss points and weights
C     and loops around all Gauss points.

      CALL GPOINTSW10(XP,XW)
      NGPT=8
      SMS=0.0D0
      SMT=0.0D0

      DO NN=1,NGPT

        RII=XP(1,NN)
        SII=XP(2,NN)
        TII=XP(3,NN)
        ALF=XW(NN)

C       Computes B matrix.
ccccccccc ESTA ES LA RUTINA CON EL PROBLEMA.
C       CALL STDM3D8(JELEM,NNODE,NDOFEL,NTENS,COORDS,MCRD,B,DDET,RII,SII,
        CALL STDM3D8(NNODE,NDOFEL,NTENS,COORDS,MCRD,B,DDET,RII,SII,
     1               TII,XBAR)

c	    CALL STDM9(NNODE,NDOFEL,NTENS,COORDS,MCRD,B,DDET,RII,SII,
c    1             XBAR)

        SMT=SMT+RO*DDET*XBAR*ALF

C       Assembles stiffness and mass matrix.

        CALL MMULT(DDSDDE,NTENS,NTENS,B,NTENS,NDOFEL,AUX1)
        CALL MTRAN(B,NTENS,NDOFEL,BT)
        CALL MMULT(BT,NDOFEL,NTENS,AUX1,NTENS,NDOFEL,FRST1)
C
C       CALL AMASS3D8(FRST2,NDOFEL,RII,SII)
        CALL AMASS3D8(FRST2,NDOFEL,RII,SII,TII)

C       Considers Gauss weight and Jacobian determinant
C       into partial stiffness and mass matrix.

        CALL SMULT(FRST1,NDOFEL,NDOFEL,ALF*DDET*XBAR)
        CALL SMULT(FRST2,NDOFEL,NDOFEL,ALF*DDET*XBAR*RO)

C       Updates stiffness and mass matrix.

        CALL UPDMAT(AMATRX,NDOFEL,NDOFEL,FRST1)
        CALL UPDMAT(AMASS ,NDOFEL,NDOFEL,FRST2)

      END DO

C     Computes elastic,inertial and damping nodal forces.

      CALL EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,
     1                    AALFA,ABETA,U,V,A,DT,SMT,RHS,JELEM)

      RETURN

      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,AALFA,ABETA,U,C
C    1                      DU,UP,DT,SMT,RHS)                          C
C                                                                      C
C   Diagonalizes the mass matrix and computes the strain and inertial  C
C   force contributions                                                C
C                                                                      C
C   AMASS(,)     :Mass matrix output in diagaonal form                 C
C   AMATRX(,)    :Stiffness matrix                                     C
C   ANMASS()     :Diagonalized mass matrix in vector form              C
C   CDMAT(,)     :Rayleigh damping matrix                              C
C   NDOFEL       :Nodal D.O.F per element                              C
C   AALFA        :Rayleigh damping parameter                           C
C   ABETA        :Rayleigh damping parameter                           C
C   U()          :Nodal displacements vector at the current increment  C
C   DU()         :Nodal incremental displacements vector               C
C   UM()         :Nodal displacements vector at the previous increment C
C   DT           :Time increment                                       C
C   SMT          :Total element mass                                   C
C   RHS()        :Right Hand Side vector                               C
C   JELEM        :Element identifier                                   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE EFFECTIVEFORCE(AMASS,AMATRX,ANMASS,NDOFEL,AALFA,
     1                          ABETA,U,V,A,DT,SMT,RHS,JELEM)

      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION AMATRX(NDOFEL,NDOFEL),AMASS(NDOFEL,NDOFEL),
     1          ANMASS(NDOFEL),U(NDOFEL),V(NDOFEL),A(NDOFEL),
     2          RSC(NDOFEL),RDC(NDOFEL),RIC(NDOFEL),RIM(NDOFEL),
     3          RSM(NDOFEL),RHS(NDOFEL)

C     Diagonalizes the mass matrix

      SMS=0.0D0
      DO I=1,NDOFEL
        SMS=SMS+AMASS(I,I)
      END DO
      DO I=1,NDOFEL
        ANMASS(I)=2.0*(AMASS(I,I)/SMS)*SMT
      END DO
C
      AMASS=0.0D0
      DO I=1,NDOFEL
        AMASS(I,I)=ANMASS(I)
      END DO

C     Forms the effective forces vector

      CALL MAVEC(AMATRX,NDOFEL,NDOFEL,U,RSC)    ! Elastic contribution
      CALL SMULTV(RSC,NDOFEL,-1.D0)
      CALL MAVEC(AMASS,NDOFEL,NDOFEL,A, RIC)    ! Inertial contribution
      CALL SMULTV(RIC,NDOFEL,-1.D0)

cc      CALL MAVEC(AMATRX,NDOFEL,NDOFEL,V,RDC)  ! Damping contribution
cc      CALL SMULTV(RDC,NDOFEL,-ABETA)
cc      CALL MAVEC(AMASS,NDOFEL,NDOFEL,V, RIM)
cc      CALL SMULTV(RIM,NDOFEL,-AALFA)

C     Puts everything into the RHS vector.

      CALL UPDVEC(RHS,NDOFEL,RSC)
      CALL UPDVEC(RHS,NDOFEL,RIC)
cc      CALL UPDVEC(RHS,NDOFEL,RDC)
cc      CALL UPDVEC(RHS,NDOFEL,RIM)

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE LINJAC(IDF,COORDS,MCRD,NNODE,IELCON,II,IM,IJ,DETJAC,    C
C    1                CT,ST)                                           C
C                                                                      C
C   Using the local connectivity array IELCON for an 8-noded 2D        C
C   element and the face ID computes the Jacobian for a 3-noded        C
C   1D-element and the outward normal vector.                          C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE LINJAC(IDF,COORDS,MCRD,NNODE,IELCON,II,IM,IJ,DETJAC,
     1                  XM,ZM,CT,ST)

      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION COORDS(MCRD,NNODE),IELCON(4,3)

      CALL LOCIEL(IELCON)

      II=IELCON(IDF,1)
      IM=IELCON(IDF,3)
      IJ=IELCON(IDF,2)

      XI=COORDS(1,II)
      ZI=COORDS(2,II)
      XM=COORDS(1,IM)
      ZM=COORDS(2,IM)
      XJ=COORDS(1,IJ)
      ZJ=COORDS(2,IJ)
      DXS=(XI-XJ)**2
      DZS=(ZI-ZJ)**2
      DL=DSQRT(DXS+DZS)
      DETJAC=DL/2.0

      CALL NORMVCTINT(XI,XJ,ZI,ZJ,CT,ST)

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE ISOCART8(RR,COORDS,NDOFDIM,XX)                          C
C                                                                      C
C   RECOVERS THE PHYSICAL CARTESIAN COORDINATES OF A POINT             C
C   THAT BELONGS TO AN 8-NODED ELEMENT GIVEN ITS NATURAL COORDINATES   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE ISOCART8(RR,COORDS,MCRD,NNODES,XX)

      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION RR(2),XX(2),COORDS(MCRD,NNODES),SN(8)

      R=RR(1)
      S=RR(2)
      CALL SHAPEF8(SN,R,S)

      XXX=0.0D0
      ZZZ=0.0D0
      DO JJ=1,8
        XXX=XXX+SN(JJ)*COORDS(1,JJ)
        ZZZ=ZZZ+SN(JJ)*COORDS(2,JJ)
      END DO

      XX(1)=XXX
      XX(2)=ZZZ

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE LOCASSEM(RHS,FF,II,IM,IJ)                               C
C                                                                      C
C   Assembles the surface force vector FF() into the elemental force   C
C   vector RHS() according to II,IM,IJ.                                C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE LOCASSEM(RHS,FF,II,IM,IJ)

      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION RHS(16),FF(6)

      RHS(2*II-1)=FF(1)
      RHS(2*II  )=FF(2)
      RHS(2*IM-1)=FF(5)
      RHS(2*IM  )=FF(6)
      RHS(2*IJ-1)=FF(3)
      RHS(2*IJ  )=FF(4)

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE RSETUP(IDF,RR)                                          C
C                                                                      C
C   According to IDF sets up the natural coordinates variation.        C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE RSETUP(IDF,RR,IIF)

      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION RR(2)

      SELECT CASE (IDF)
        CASE (1)
        RR(2)=-1.0
        IIF=1

        CASE (2)
        RR(1)= 1.0D0
        IIF=2

        CASE (3)
        RR(2)= 1.0D0
        IIF=1

        CASE (4)
        RR(1)=-1.0D0
        IIF=2

      END SELECT

      RETURN

      END

ext: for


text:
submodule (vtk_element) vtk_element_procedures
    implicit none
    !! author: Ian Porter
    !! date: 05/06/2019
    !!
    !! this is the basic file for a serial vtk file
    !!

contains

    module procedure vtk_element_setup
        use file_utility, only : is_little_endian
        implicit none
        !! author: Ian Porter
        !! date: 05/06/2019
        !!
        !! this writes the header for a vtk file
        !!
        !! example:
        !! <VTKFile type=”ImageData” version=”0.1” byte_order=”LittleEndian”>
        character(len=*), parameter   :: name = 'VTKFile'
        character(len=:), allocatable :: string
        character(len=:), allocatable :: type_string
        character(len=*), parameter   :: version_string = ' version="' // def_version // '"'
        character(len=:), allocatable :: byte_order_string
        character(len=:), allocatable :: compression_string

        if (allocated(me%type)) then
            allocate(type_string,source=' type="' // me%type // '"')
        else
            error stop "error. can't create vtk file without a known type. terminated in vtk_element_setup"
        end if
        if (is_little_endian()) then
            allocate(byte_order_string,source=' byte_order="LittleEndian"')
        else
            allocate(byte_order_string,source=' byte_order="BigEndian"')
        end if
        if (allocated(me%compression)) then
            allocate(compression_string,source=' compression="' // me%compression // '"')
        else
            allocate(compression_string,source='')
        end if

        allocate(string, source=type_string // version_string // byte_order_string // compression_string)

        call me%setup(name=name,string=string,offset=4)

    end procedure vtk_element_setup

    module procedure initialize
        implicit none
        !! author: Ian Porter
        !! date: 05/07/2019
        !!
        !! this is an external interface to allow the user or other routines to set the internal variables
        !!

        if (present(type))           allocate(me%type,source=type)
        if (present(compression))    allocate(me%compression,source=compression)
        if (present(file_extension)) allocate(me%file_extension,source=file_extension)

        call me%vtk_element_setup()

    end procedure initialize

    module procedure finalize
        implicit none
        !! author: Ian Porter
        !! date: 05/07/2019
        !!
        !!
        !!
        if (allocated(me%vtk_element)) then
            call me%vtk_element%finalize()
            call me%add(me%vtk_element)
        end if

    end procedure finalize

    module procedure gcc_bug_workaround_deallocate_vtk_element_single
        implicit none
        !! gcc work-around for deallocating a multi-dimension derived type w/ allocatable character strings

        if (allocated(foo%type))           deallocate(foo%type)
        if (allocated(foo%version))        deallocate(foo%version)
        if (allocated(foo%compression))    deallocate(foo%compression)
        if (allocated(foo%file_extension)) deallocate(foo%file_extension)
        if (allocated(foo%filename))       deallocate(foo%filename)
        !call foo%piece%deallocate_piece_dt()
        if (allocated(foo%vtk_element)) then
            call foo%vtk_element%me_deallocate()
            call foo%vtk_element%deallocate()
        end if
        call foo%deallocate()

    end procedure gcc_bug_workaround_deallocate_vtk_element_single

end submodule vtk_element_procedures

max_stars_repo_name: porteri/vtkmofo

max_issues_repo_name: porteri/vtkmofo

max_forks_repo_name: nrc-fuels/vtkmofo