text:
Blue Oak Energys mission is to transform the ways we generate electricity by making sustainable solar energy widely accessible. Blue Oak Energy has engineered and constructed solar electric systems on over 600 distinct locations across more than 30 states in the US. The work they have completed has added more than 500MW of clean renewable energy to the utility grid which will produce energy for 30 years or more. Their success is a direct result of our dedicated and talented technical team who work solely on commercial and utility solar facilities.

Recent projects include a comprehensive roofmounted PV and solar carport design package for the Google Campus headquarters in Mountain View.

They were formerly located at 4614 Second Street. The move allowed them to expand in size.

 see Solar Power

 Media
  http://www.davisenterprise.com/business/solarcompanymovesitsdavisheadquarters/ Solar company moves its Davis headquarters
  http://www.businesswire.com/news/home/20130402005547/en/BlueOakEnergyBuilds4.2MWUtility#.Uym1E_ldVgk
  http://www.blueoakenergy.com/documents/presskit/pressreleases/Inc5000_Press_Release_for_Website_2013.pdf


text:
subroutine depfil_double(lundia    ,error     ,fildep    ,fmttmp    ,array     , &
                       & nfld      ,ifld      ,dims      )
!----- LGPL --------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This library is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License as published by the Free Software Foundation version 2.1.
!
!  This library is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with this library; if not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: depfil_double.f90 5915 2016-03-02 10:39:57Z jagers $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/utils_lgpl/deltares_common/packages/deltares_common_mpi/src/depfil_double.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Reads the depth values from the attribute file
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use grid_dimens_module
    use message_module
    use system_utils, only: exifil
    use quickin_dep_file
    use dfparall
    !
    implicit none
    !
    ! The following list of pointer parameters is used to point inside the dims structure
    !
    integer, pointer :: mfg
    integer, pointer :: mlg
    integer, pointer :: nfg
    integer, pointer :: nlg
    integer, pointer :: mmaxgl
    integer, pointer :: nmaxgl
!
! Global variables
!
    type(griddimtype), target                                          , intent(in)  :: dims   !  grid dimensions
    integer                                                            , intent(in)  :: ifld   !  index of field to be read
    integer                                                                          :: lundia !  unit number for diagnostic file
    integer                                                            , intent(in)  :: nfld   !  number of fields
    logical                                                            , intent(out) :: error  !  Flag=TRUE if an error is encountered
    real(hp), dimension(nfld, dims%nlb:dims%nub, dims%mlb:dims%mub)    , intent(out) :: array  !  data array to fill
    character(*)                                                                     :: fildep !  Name of the relevant file
    character(11)                                                      , intent(in)  :: fmttmp !  Format switch for the attribute file
!
! Local variables
!
    integer                               :: iocond ! Help variable for iostat condition
    integer                               :: luntmp ! Unit number for attribute file
    integer                               :: m
    integer                               :: n
    integer                 , external    :: newunit
    real(hp), dimension(:,:), allocatable :: dtmp   ! Temporary array containing dp of entire domain
    character(300)                        :: errmsg ! Character string containing the errormessage to be written to file. The message depends on the error.
    !
    integer                               :: i
    integer                               :: nmaxio
    integer                               :: mmaxio
    integer , dimension(:,:), allocatable :: cnt
    real(hp), dimension(:,:), allocatable :: dtmpio
!
!! executable statements -------------------------------------------------------
!
    mfg    => dims%mfg
    mlg    => dims%mlg
    nfg    => dims%nfg
    nlg    => dims%nlg
    mmaxgl => dims%mmaxgl
    nmaxgl => dims%nmaxgl
    !
    error = .false.
    !
    ! Test file existence, if so read
    !
    if (exifil(fildep, lundia)) then
       !
       ! File exists
       !
       ! allocate temporary array to store data of entire domain read from file
       !
       ! NOTE: nmaxus and mmax equal nmaxgl and mmaxgl, respectively (for entire domain)
       !       in case of parallel runs. Moreover, array is associated with subdomain and
       !       therefore, data for entire domain is stored in temporary array dtmp
       !
       allocate (dtmp(nmaxgl,mmaxgl))
       !
       ! the master opens and reads the depth file
       !
       if ( inode /= master ) goto 10
       !
       if (associated(dims%aggrtable)) then
          mmaxio = size(dims%aggrtable,1)
          nmaxio = size(dims%aggrtable,2)
          allocate (dtmpio(nmaxio,mmaxio),cnt(nmaxgl,mmaxgl))
          !
          call read_qndep(lundia    ,error     ,fildep    ,fmttmp    ,dtmpio    , &
                        & nmaxio    ,mmaxio    ,iocond    )
          cnt = 0
          dtmp = 0.0_hp
          do m = 1, mmaxio
             do n = 1, nmaxio
                i = dims%aggrtable(n,m)
                if (i>0) then
                   dtmp(i,1) = dtmp(i,1) + dtmpio(n,m)
                   cnt(i,1) = cnt(i,1) + 1
                endif
             enddo
          enddo
          do i = 1, nmaxgl
             if (cnt(i,1)>1) then
                dtmp(i,1) = dtmp(i,1)/cnt(i,1)
             endif
          enddo
          !
          deallocate(dtmpio,cnt)
       else
          call read_qndep(lundia    ,error     ,fildep    ,fmttmp    ,dtmp      , &
                        & nmaxgl    ,mmaxgl    ,iocond    )
       endif
       !
 10    continue
       !
       ! check whether something went wrong with reading file
       !
       call dfbroadc ( iocond, 1, dfint, error, errmsg )
       if ( error ) then
          !
          ! show error message when not yet done so
          !
          if (errmsg /= ' ') then
             call write_error(errmsg, unit=lundia)
          endif
       elseif (iocond /= 0) then
          if (iocond < 0) then
             errmsg = PREMATURE_EOF // trim(fildep)
             call write_error(errmsg, unit=lundia)
          else
             errmsg = FILE_READ_ERROR // trim(fildep)
             call write_error(errmsg, unit=lundia)
          endif
          error = .true.
       else
          !
          ! scatter array dtmp to all nodes
          !
          call dfbroadc ( dtmp, nmaxgl*mmaxgl, dfdble, error, errmsg )
          if ( error ) then
             call write_error(errmsg, unit=lundia)
          else
             !
             ! put copies of parts of dtmp for each subdomain
             !
             do m = mfg, mlg
                do n = nfg, nlg
                   array(ifld,n-nfg+1,m-mfg+1) = dtmp(n,m)
                enddo
             enddo
          endif
       endif
       !
       deallocate(dtmp)
    else
       !
       ! File does not exist
       ! Exifil has produced a nice error message
       !
       error = .true.
    endif
end subroutine depfil_double

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      program symm
      implicit none

      integer nosp,i,j,k,n,nsea,toteq,nobsea,ll,ijk,is,kk,js,l,lp
      real(8) diff,dumx,dumy,dumz,anorm,df1,df2,df3,df4,norm1,norm2
      real(8) anx2,any2,anz2,dp
      real(8),dimension(:),allocatable:: sai
      real(8),dimension(:,:),allocatable:: d
      real(8),dimension(:,:),allocatable:: x,xf
      real(8),dimension(3)   :: dx,vd1,vd2
      real(8),dimension(3,3) :: sig
      integer,dimension(:),allocatable:: ian,nn,icode,nn2
      integer,dimension(:,:),allocatable:: isea,isea2

      read(5,*) n
      allocate(sai(n))
      allocate(ian(n),nn(n),icode(n),nn2(n))
      allocate(d(n,n))
      allocate(isea(n,n),isea2(n,n))
      allocate(x(3,n),xf(3,n))

      nosp=0
      nobsea=0
      do i=1,n
         nn(i)=0
         read(5,*) ian(i),x(1,i),x(2,i),x(3,i)
      enddo
      nsea=0
c construct the distance matrix
      do i=1,n
         sai(i)=0
         do j=1,n
            dumx=x(1,i)-x(1,j)
            dumy=x(2,i)-x(2,j)
            dumz=x(3,i)-x(3,j)
            d(i,j)=dsqrt(dumx*dumx+dumy*dumy+dumz*dumz)
            sai(i)=sai(i)+ian(i)*d(i,j)**3
         enddo
      enddo

      do i=1,n
         icode(i)=0
         do k=1,n
            diff=abs(sai(i)-sai(k))
c            print*, i,k,sai(i),sai(k),diff
            if(diff<=0.15d0) then
               nn(i)=nn(i)+1
               isea(i,nn(i))=k
               icode(i)=icode(i)+k**k
            endif
         enddo
         if(nn(i)>1) then
            lp=1
            do l=1,i-1
               if(icode(i)==icode(l)) lp=0
            enddo
            if(lp==1)  then
               nobsea=nobsea+1
               nn2(nobsea)=nn(i)
               do ll=1,nn2(nobsea)
                  isea2(nobsea,ll)=isea(i,ll)
c                  print*, nobsea,ll,isea2(nobsea,ll)
               enddo
            endif
         endif
      enddo
      do i=1,nobsea
c         print*, "i",i,nobsea
         do j=1,nn2(i)
            do k=j+1,nn2(i)
c all possible j-k combinations
c               print*, i,j,k,isea2(i,j),isea(i,k)
               dx(1)=x(1,isea2(i,j))-x(1,isea2(i,k))
               dx(2)=x(2,isea2(i,j))-x(2,isea2(i,k))
               dx(3)=x(3,isea2(i,j))-x(3,isea2(i,k))
               anx2=dx(1)*dx(1)
               any2=dx(2)*dx(2)
               anz2=dx(3)*dx(3)
               anorm=dsqrt( anx2+any2+anz2 )
               dx(1)=dx(1)/anorm
               dx(2)=dx(2)/anorm
               dx(3)=dx(3)/anorm
               do is=1,3
                  do js=is,3
                     if(is==js) then
                        sig(is,js)=1-2*dx(is)*dx(js)
                     else
                        sig(is,js)=-2*dx(is)*dx(js)
                        sig(js,is)=sig(is,js)
                     endif
                  enddo
               enddo
               xf=matmul(sig,x)
               toteq=0
               do ijk=1,n
                  do kk=1,n
                     df1=ian(ijk)-ian(kk)
                     df2=xf(1,ijk)-x(1,kk)
                     df3=xf(2,ijk)-x(2,kk)
                     df4=xf(3,ijk)-x(3,kk)
                     diff=df1*df1+df2*df2+df3*df3+df4*df4
                     if(diff<=0.0001d0) toteq=toteq+1
                  enddo
                  if(toteq==n) then
                      nosp=nosp+1
                      print*, "Normal direction of symm. plane",
     &                isea2(i,j),isea2(i,k)
                  endif
               enddo
            enddo
         enddo
      enddo
      print*, "N of symmetry planes",nosp
      if(nosp>0)  print*, "Symmetry: CS"
      if(nosp==0) print*, "Symmetry: C1"
1     format('H',3(f20.10))
6     format('C',3(f20.10))
      end

ext: f

max_stars_repo_name: emartineznunez/AutoMeKin

max_issues_repo_name: emartineznunez/AutoMeKin

max_forks_repo_name: emartineznunez/AutoMeKin


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!
!     Determining the coefficients of the dual shape functions on
!     the slave surface
!     see phd thesis Sitzmann Chapter 3.3. (p.34) and Chapter 4.3.
!
!
!     [in]     islavact  (i) indicates, if slave node i is active (=-3 no-slave-node, =-2 no-LM-node, =-1 no-gap-node, =0 inactive node, =1 sticky node, =2 slipping/active node)
!     [in]     islavsurf  islavsurf(1,i) slaveface i islavsurf(2,i) pointer into imastsurf and pmastsurf
!     [in]     itiefac    pointer into field islavsurf: (1,i) beginning slave_i (2,i) end of slave_i
!     [in]     islavnode  fields containing nodes of slace surfaces
!     [in]     nslavnode  (i) for contraint i pointer into field islavnode
!     [in]     pslavsurf  integration points and weights on slave side
!     [in] pslavdual (:,i)coefficients \f$ \alpha_{ij}\f$, \f$ 1,j=1,..8\f$ for dual shape functions for face i
!     [in] pslavdualpg (:,i)coefficients \f$ \alpha_{ij}\f$, \f$ 1,j=1,..8\f$ for Petrov-Galerkin shape functions for face i
!     [in]  iflagdualquad   flag indicating what mortar contact is used (=1 quad-lin, =2 quad-quad, =3 PG quad-lin, =4 PG quad-quad)
!
      subroutine gendualcoeffs(tieset,ntie,ipkon,kon,lakon,co,vold,
     &     islavact,islavsurf,itiefac,islavnode,nslavnode,
     &     mi,pslavsurf,pslavdual,pslavdualpg,iflagdualquad)
!
!     Determining the coefficients of the dual shape functions on
!     the slave surface
!
!     Author: Sitzmann,Saskia ;
!
      implicit none
!
      logical checkbiorthogonality, checknorm
!
      character*8 lakon(*)
      character*81 tieset(3,*),slavset
!
      integer ntie,ipkon(*),kon(*),iflag,kneigh,i,ii,j,k,l,
     &     ipos,islavact(*),nope,islavsurf(2,*),islavnode(*),
     &     nslavnode(ntie+1),itiefac(2,*),ifaces,nelems,jfaces,mi(*),
     &     mint2d,m,nopes,konl(20),id,nopes2,ipiv(8),info,ipnt,ifac,
     &     getlocno,lnode(2,8),nnogap,n1,n2,iscontr(8*8),
     &     imcontr(8*8),jj,locm,locs,nodesf,nodem,ifs,
     &     ifm,ns,indexf,idummy,modf,iflagdualquad
!
      real*8 co(3,*),vold(0:mi(2),*),ets,xis,weight,xl2s(3,8),xsj2(3),
     &     shp2(7,8),shp2s(7,8),dx,help,xs2(3,7),
     &     pslavdual(64,*),diag_els(8),m_els(36),contribution,work(8),
     &     contr(8*8),xs2m(3,7),xsj2m(3),shp2m(7,8),etm,xim,
     &     pslavsurf(3,*),pslavdualpg(64,*)
!
!
!
      data iflag /2/
!
      checkbiorthogonality=.false.
      checknorm=.false.
!
!     loop over all ties
!
      do i=1,ntie
        if(tieset(1,i)(81:81).ne.'C') cycle
        kneigh=1
        slavset=tieset(2,i)
        ipos=index(slavset,' ')
        if(slavset(ipos-1:ipos-1).eq.'S') cycle
!
!     ntri: number of triangles in the triangulation of the master
!     surface corresponding to tie i
!
        do l=itiefac(1,i),itiefac(2,i)
          if(checkbiorthogonality)then
            write(*,*)'face',l
          endif
          ifaces=islavsurf(1,l)
          nelems=int(ifaces/10)
          jfaces=ifaces-nelems*10
          call getnumberofnodes(nelems,jfaces,lakon,nope,
     &         nopes,idummy)
!
!     initialization for Dualshape Coefficient matrix
!
          do j=1,8
            diag_els(j)=0.0
            do k=1,j
              m_els(k+((j-1)*j/2))=0.0
            enddo
          enddo
!
!     actual position of the nodes belonging to the
!     slave surface
!
          do j=1,nope
            konl(j)=kon(ipkon(nelems)+j)
          enddo
          nnogap=0
!
!     check for nogap or nolm Nodes
!     here redistribution of Lagrange Multiplier contribution
!     has to be applied
!     see Phd thesis Sitzmann Chapter 4.3.
!
          do m=1,nopes
            ifac=getlocno(m,jfaces,nope)
            lnode(1,m)=konl(ifac)
            call nident(islavnode(nslavnode(i)+1),
     &           konl(ifac),(nslavnode(i+1)-nslavnode(i)),id)
            if(islavnode(nslavnode(i)+id).eq.konl(ifac)) then
              lnode(2,m)=islavact(nslavnode(i)+id)
!
!     for quad-lin mortar method only corner nodes carry
!     Lagrange Multiplier contribution
!
              if((nopes.eq.8).and.
     &             ((iflagdualquad.eq.1).or.(iflagdualquad.eq.3)))
     &             then
                if((lnode(2,m).lt.0).and.(m.le.4))then
                  nnogap=nnogap+1
                endif
                if(m.gt.4) then
                  islavact(nslavnode(i)+id)=-2
                  lnode(2,m)=islavact(nslavnode(i)+id)
                endif
              else if((nopes.eq.6).and.
     &               ((iflagdualquad.eq.1).or.(iflagdualquad.eq.3)))then
                if((lnode(2,m).lt.0).and.(m.le.3)) then
                  nnogap=nnogap+1
                endif
                if(m.gt.3) then
                  islavact(nslavnode(i)+id)=-2
                  lnode(2,m)=islavact(nslavnode(i)+id)
                endif
              else
                if(lnode(2,m).lt.0) nnogap=nnogap+1
              endif
              if(checknorm)then
                write(*,*) 'node',lnode(1,m),lnode(2,m)
              endif
            else
              write(*,*)'createbd: node',konl(ifac)
              write(*,*)'was not catalogued properly in islavnode'
              call exit(201)
            endif

            do j=1,3
              xl2s(j,m)=co(j,konl(ifac))+
     &             vold(j,konl(ifac))
            enddo
          enddo
          if(((iflagdualquad.eq.1).or.(iflagdualquad.eq.3)))then
            if(nopes.eq.6)then
              nopes2=3
            else if(nopes.eq.8)then
              nopes2=4
            else
              nopes2=nopes
            endif
          else
            nopes2=nopes
          endif
!
          mint2d=islavsurf(2,l+1)-islavsurf(2,l)
          if(mint2d.eq.0) cycle
          indexf=islavsurf(2,l)
!
!     Calculate the Mass matrix for compilation of the dualshapefunction
!     pslavdual(16,*)
!     create Z, W (see Phd thesis Sitzmann equation (3.34))
!
!     loop over integration points
!
          do m=1,mint2d
            xis=pslavsurf(1,indexf+m)
            ets=pslavsurf(2,indexf+m)
            weight=pslavsurf(3,indexf+m)
            ns=l
            iflag=2
!
!     construct dual basis functions for transformed standart
!     basis functions
!
            if(nopes.eq.8) then
              if((iflagdualquad.eq.2).or.(iflagdualquad.eq.4))then
                call shape8qtilde(xis,ets,xl2s,xsj2,xs2,shp2,iflag)
              else
                call shape8qtilde_lin(xis,ets,xl2s,xsj2,xs2,
     &               shp2,iflag)
              endif
            elseif(nopes.eq.4) then
              call shape4q(xis,ets,xl2s,xsj2,xs2,shp2,iflag)
            elseif(nopes.eq.6) then
              if((iflagdualquad.eq.2).or.(iflagdualquad.eq.4))then
                call shape6tritilde(xis,ets,xl2s,xsj2,xs2,shp2,
     &               iflag)
              else
                call shape6tritilde_lin(xis,ets,xl2s,xsj2,xs2,
     &               shp2,iflag)
              endif
            else
              call shape3tri(xis,ets,xl2s,xsj2,xs2,shp2,iflag)
            endif
            dx=dsqrt(xsj2(1)**2+xsj2(2)**2+
     &           xsj2(3)**2)
            ipnt=0
!
!     column
!
            do j=1,nopes2
              diag_els(j)=diag_els(j)+shp2(4,j)*weight
     &             *dx
              do k=1,j
                m_els(k+((j-1)*j/2))=m_els(k+((j-1)*j/2))
     &               +shp2(4,k)*shp2(4,j)
     &               *weight
     &               *dx
              enddo
            enddo
          enddo
!
!     compute inverse of me_ls
!     factorisation
!
          if(checknorm .or.checkbiorthogonality)then
            write(*,*)'diag_els',l
            write(*,105)(diag_els(j),j=1,nopes2)
            write(*,*)'m_els'
            do j=1,nopes2
              write(*,105) (m_els(k+((j-1)*j/2)),k=1,j)
            enddo
          endif
          call dsptrf('U',nopes2,m_els,ipiv,info)
!
!     inverse
!
          call dsptri('U',nopes2,m_els,ipiv,work,info)
!
!     stack of pslavdual multiplication with diag_els
!     A=W*Z^-1
!
          do k=1,64
            pslavdual(k,l)=0.0
            if(iflagdualquad.gt.2)then
              pslavdualpg(k,l)=0.0
            endif
          enddo
          do k=1,nopes2
            if(iflagdualquad.gt.2)then
              pslavdualpg((k-1)*8+k,l)=1.0
            endif
            do j=1,nopes2
              if(k.le.j)then
                pslavdual((k-1)*8+j,l)=diag_els(k)
     &               *m_els(k+((j-1)*j/2))
              else
                pslavdual((k-1)*8+j,l)=diag_els(k)
     &               *m_els(j+((k-1)*k/2))
              endif
            enddo
          enddo
!
          if(checknorm)then
            write(*,*)'pslavdual1,face',l
            do ii=1,nopes2
              help=0.0
              do j=1,nopes2
                help=help+pslavdual((ii-1)*8+j,l)
              enddo
              write(*,105) (pslavdual((ii-1)*8+j,l),j=1,nopes2)
            enddo
          endif
!
!     remove Lagrange Multiplier contribution for nogap and noLM nodes
!     phd thesis Sitzmann Chapter 4.3.
!
          if((nnogap.gt.0).and.(nopes2.eq.4)) then
            if(nnogap.eq.1) then
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) exit
              enddo
              n1=modf(nopes2,ii-1)
              n2=modf(nopes2,ii+1)
              do jj=1,nopes2
                pslavdual((n1-1)*8+jj,l)=pslavdual((n1-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                pslavdual((n2-1)*8+jj,l)=pslavdual((n2-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                if(iflagdualquad.gt.2)then
                  pslavdualpg((n1-1)*8+jj,l)=
     &                 pslavdualpg((n1-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                  pslavdualpg((n2-1)*8+jj,l)=
     &                 pslavdualpg((n2-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                endif
              enddo
              do jj=1,nopes2
                pslavdual((ii-1)*8+jj,l)=0.0
                if(iflagdualquad.gt.2)then
                  pslavdualpg((ii-1)*8+jj,l)=0.0
                endif
              enddo
            elseif(nnogap.eq.2)then
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) then
                  n1=modf(nopes2,ii-1)
                  n2=modf(nopes2,ii+1)
                  if(lnode(2,n1).lt.0) n1=n2
                  if(lnode(2,n2).lt.0) n2=n1
                  do jj=1,nopes2
                    pslavdual((n1-1)*8+jj,l)=
     &                   pslavdual((n1-1)*8+jj,l)
     &                   +0.5*pslavdual((ii-1)*8+jj,l)
                    pslavdual((n2-1)*8+jj,l)=
     &                   pslavdual((n2-1)*8+jj,l)
     &                   +0.5*pslavdual((ii-1)*8+jj,l)
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((n1-1)*8+jj,l)=
     &                     pslavdualpg((n1-1)*8+jj,l)
     &                     +0.5*pslavdualpg((ii-1)*8+jj,l)
                      pslavdualpg((n2-1)*8+jj,l)=
     &                     pslavdualpg((n2-1)*8+jj,l)
     &                     +0.5*pslavdualpg((ii-1)*8+jj,l)
                    endif
                  enddo
                  do jj=1,nopes2
                    pslavdual((ii-1)*8+jj,l)=0.0
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((ii-1)*8+jj,l)=0.0
                    endif
                  enddo
                endif
              enddo
            elseif(nnogap.eq.3)then
              do ii=1,nopes2
                if(lnode(2,ii).gt.-1) n1=ii
              enddo
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) then
                  do jj=1,nopes2
                    pslavdual((n1-1)*8+jj,l)=
     &                   pslavdual((n1-1)*8+jj,l)
     &                   +pslavdual((ii-1)*8+jj,l)
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((n1-1)*8+jj,l)=
     &                     pslavdualpg((n1-1)*8+jj,l)
     &                     +pslavdualpg((ii-1)*8+jj,l)
                    endif
                  enddo
                  do jj=1,nopes2
                    pslavdual((ii-1)*8+jj,l)=0.0
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((ii-1)*8+jj,l)=0.0
                    endif
                  enddo
                endif
              enddo
            endif
          endif
          if((nnogap.gt.0).and.(nopes2.eq.3)) then
            if(nnogap.eq.1) then
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) exit
              enddo
              n1=modf(nopes2,ii-1)
              n2=modf(nopes2,ii+1)
              do jj=1,nopes2
                pslavdual((n1-1)*8+jj,l)=pslavdual((n1-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                pslavdual((n2-1)*8+jj,l)=pslavdual((n2-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                if(iflagdualquad.gt.2)then
                  pslavdualpg((n1-1)*8+jj,l)=
     &                 pslavdualpg((n1-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                  pslavdualpg((n2-1)*8+jj,l)=
     &                 pslavdualpg((n2-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                endif
              enddo
              do jj=1,nopes2
                pslavdual((ii-1)*8+jj,l)=0.0
                if(iflagdualquad.gt.2)then
                  pslavdualpg((ii-1)*8+jj,l)=0.0
                endif
              enddo
            elseif(nnogap.eq.2)then
              do ii=1,nopes2
                if(lnode(2,ii).gt.-1) n1=ii
              enddo
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) then
                  do jj=1,nopes2
                    pslavdual((n1-1)*8+jj,l)=
     &                   pslavdual((n1-1)*8+jj,l)
     &                   +pslavdual((ii-1)*8+jj,l)
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((n1-1)*8+jj,l)=
     &                     pslavdualpg((n1-1)*8+jj,l)
     &                     +pslavdualpg((ii-1)*8+jj,l)
                    endif
                  enddo
                  do jj=1,nopes2
                    pslavdual((ii-1)*8+jj,l)=0.0
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((ii-1)*8+jj,l)=0.0
                    endif
                  enddo
                endif
              enddo
            endif
          elseif((nnogap.gt.0).and.(nopes2.eq.8))then
            if(nnogap.eq.1)then
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) exit
              enddo
              if(ii.le.4)then
                n1=modf(4,ii-1)+4
                n2=ii+4
              else
                n1=ii-4
                n2=modf(4,ii-4+1)
              endif
              do jj=1,nopes2
                pslavdual((n1-1)*8+jj,l)=pslavdual((n1-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                pslavdual((n2-1)*8+jj,l)=pslavdual((n2-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                if(iflagdualquad.gt.2)then
                  pslavdualpg((n1-1)*8+jj,l)=
     &                 pslavdualpg((n1-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                  pslavdualpg((n2-1)*8+jj,l)=
     &                 pslavdualpg((n2-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                endif
              enddo
              do jj=1,nopes2
                pslavdual((ii-1)*8+jj,l)=0.0
                if(iflagdualquad.gt.2)then
                  pslavdualpg((ii-1)*8+jj,l)=0.0
                endif
              enddo
            elseif((nnogap.gt.1).and.(nnogap.lt.7))then
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) then
                  if(ii.le.4)then
                    n1=modf(4,ii-1)+4
                    n2=ii+4
                  else
                    n1=ii-4
                    n2=modf(4,ii-4+1)
                  endif
 10               if((lnode(2,n1).lt.0).and.(lnode(2,n2).lt.0))then
                    if(n1.le.4)then
                      n1=modf(4,n1-1)+4
                    else
                      n1=n1-4
                    endif
                    if(n2.le.4)then
                      n2=n2+4
                    else
                      n2=modf(4,n2-4+1)
                    endif
                    go to 10
                  else
                    if(lnode(2,n1).lt.0) n1=n2
                    if(lnode(2,n2).lt.0) n2=n1
                    do jj=1,nopes2
                      pslavdual((n1-1)*8+jj,l)=
     &                     pslavdual((n1-1)*8+jj,l)
     &                     +0.5*pslavdual((ii-1)*8+jj,l)
                      pslavdual((n2-1)*8+jj,l)=
     &                     pslavdual((n2-1)*8+jj,l)
     &                     +0.5*pslavdual((ii-1)*8+jj,l)
                      if(iflagdualquad.gt.2)then
                        pslavdualpg((n1-1)*8+jj,l)=
     &                       pslavdualpg((n1-1)*8+jj,l)
     &                       +0.5*pslavdualpg((ii-1)*8+jj,l)
                        pslavdualpg((n2-1)*8+jj,l)=
     &                       pslavdualpg((n2-1)*8+jj,l)
     &                       +0.5*pslavdualpg((ii-1)*8+jj,l)
                      endif
                    enddo
                  endif
                  do jj=1,nopes2
                    pslavdual((ii-1)*8+jj,l)=0.0
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((ii-1)*8+jj,l)=0.0
                    endif
                  enddo
                endif
              enddo
            elseif(nnogap.eq.7)then
              do ii=1,nopes2
                if(lnode(2,ii).gt.-1) n1=ii
              enddo
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) then
                  do jj=1,nopes2
                    pslavdual((n1-1)*8+jj,l)=
     &                   pslavdual((n1-1)*8+jj,l)
     &                   +pslavdual((ii-1)*8+jj,l)
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((n1-1)*8+jj,l)=
     &                     pslavdualpg((n1-1)*8+jj,l)
     &                     +pslavdualpg((ii-1)*8+jj,l)
                    endif
                  enddo
                  do jj=1,nopes2
                    pslavdual((ii-1)*8+jj,l)=0.0
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((ii-1)*8+jj,l)=0.0
                    endif
                  enddo
                endif
              enddo
            endif
          elseif((nnogap.gt.0).and.(nopes2.eq.6))then
            if(nnogap.eq.1)then
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) exit
              enddo
              if(ii.le.3)then
                n1=modf(3,ii-1)+3
                n2=ii+3
              else
                n1=ii-3
                n2=modf(3,ii-3+1)
              endif
              do jj=1,nopes2
                pslavdual((n1-1)*8+jj,l)=pslavdual((n1-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                pslavdual((n2-1)*8+jj,l)=pslavdual((n2-1)*8+jj,l)
     &               +0.5*pslavdual((ii-1)*8+jj,l)
                if(iflagdualquad.gt.2)then
                  pslavdualpg((n1-1)*8+jj,l)=
     &                 pslavdualpg((n1-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                  pslavdualpg((n2-1)*8+jj,l)=
     &                 pslavdualpg((n2-1)*8+jj,l)
     &                 +0.5*pslavdualpg((ii-1)*8+jj,l)
                endif
              enddo
              do jj=1,nopes2
                pslavdual((ii-1)*8+jj,l)=0.0
                if(iflagdualquad.gt.2)then
                  pslavdualpg((ii-1)*8+jj,l)=0.0
                endif
              enddo
            elseif((nnogap.gt.1).and.(nnogap.lt.5))then
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) then
                  if(ii.le.3)then
                    n1=modf(3,ii-1)+3
                    n2=ii+3
                  else
                    n1=ii-3
                    n2=modf(3,ii-3+1)
                  endif
 20               if((lnode(2,n1).lt.0).and.(lnode(2,n2).lt.0))then
                    if(n1.le.3)then
                      n1=modf(3,n1-1)+3
                    else
                      n1=n1-3
                    endif
                    if(n2.le.3)then
                      n2=n2+3
                    else
                      n2=modf(3,n2-3+1)
                    endif
                    go to 20
                  else
                    if(lnode(2,n1).lt.0) n1=n2
                    if(lnode(2,n2).lt.0) n2=n1
                    do jj=1,nopes2
                      pslavdual((n1-1)*8+jj,l)=
     &                     pslavdual((n1-1)*8+jj,l)
     &                     +0.5*pslavdual((ii-1)*8+jj,l)
                      pslavdual((n2-1)*8+jj,l)=
     &                     pslavdual((n2-1)*8+jj,l)
     &                     +0.5*pslavdual((ii-1)*8+jj,l)
                      if(iflagdualquad.gt.2)then
                        pslavdualpg((n1-1)*8+jj,l)=
     &                       pslavdualpg((n1-1)*8+jj,l)
     &                       +0.5*pslavdualpg((ii-1)*8+jj,l)
                        pslavdualpg((n2-1)*8+jj,l)=
     &                       pslavdualpg((n2-1)*8+jj,l)
     &                       +0.5*pslavdualpg((ii-1)*8+jj,l)
                      endif
                    enddo
                  endif
                  do jj=1,nopes2
                    pslavdual((ii-1)*8+jj,l)=0.0
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((ii-1)*8+jj,l)=0.0
                    endif
                  enddo
                endif
              enddo
            elseif(nnogap.eq.5)then
              do ii=1,nopes2
                if(lnode(2,ii).gt.-1) n1=ii
              enddo
              do ii=1,nopes2
                if(lnode(2,ii).lt.0) then
                  do jj=1,nopes2
                    pslavdual((n1-1)*8+jj,l)=
     &                   pslavdual((n1-1)*8+jj,l)
     &                   +pslavdual((ii-1)*8+jj,l)
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((n1-1)*8+jj,l)=
     &                     pslavdualpg((n1-1)*8+jj,l)
     &                     +pslavdualpg((ii-1)*8+jj,l)
                    endif
                  enddo
                  do jj=1,nopes2
                    pslavdual((ii-1)*8+jj,l)=0.0
                    if(iflagdualquad.gt.2)then
                      pslavdualpg((ii-1)*8+jj,l)=0.0
                    endif
                  enddo
                endif
              enddo
            endif
          endif
          if(checknorm)then
            write(*,*)'pslavdual2,face',l
            do ii=1,8
              help=0.0
              do j=1,8
                help=help+pslavdual((ii-1)*8+j,l)
              enddo
              write(*,105) (pslavdual((ii-1)*8+j,l),j=1,8)
            enddo
          endif
 105      format(8(1x,e15.4))
        enddo
!
!     FIRST SLAVE SURFACE LOOP DONE
!
        if(checkbiorthogonality)then
          do l=itiefac(1,i),itiefac(2,i)
            ifaces=islavsurf(1,l)
            nelems=int(ifaces/10)
            jfaces=ifaces-nelems*10
            call getnumberofnodes(nelems,jfaces,lakon,nope,
     &           nopes,idummy)
!
!     actual position of the nodes belonging to the
!     slave surface
!
            do j=1,nope
              konl(j)=kon(ipkon(nelems)+j)
            enddo
            do m=1,nopes
              ifac=getlocno(m,jfaces,nope)
              do j=1,3
                xl2s(j,m)=co(j,konl(ifac))+
     &               vold(j,konl(ifac))
              enddo
            enddo
            mint2d=islavsurf(2,l+1)-islavsurf(2,l)
            if(mint2d.eq.0) cycle
            indexf=islavsurf(2,l)
            do m=1,mint2d
              xis=pslavsurf(1,indexf+m)
              ets=pslavsurf(2,indexf+m)
              weight=pslavsurf(3,indexf+m)
              ns=l
              iflag=2
              if(nopes.eq.8) then
                if((iflagdualquad.eq.2).or.(iflagdualquad.eq.4))then
                  call dualshape8qtilde(xis,ets,xl2s,xsj2,xs2,
     &                 shp2s,ns,pslavdual,iflag)
                else
                  call dualshape8qtilde_lin(xis,ets,xl2s,xsj2,
     &                 xs2,shp2s,ns,pslavdual,iflag)
                endif
              elseif(nopes.eq.4) then
                call dualshape4q(xis,ets,xl2s,xsj2,xs2,shp2s,ns,
     &               pslavdual,iflag)
              elseif(nopes.eq.6) then
                if((iflagdualquad.eq.2).or.(iflagdualquad.eq.4))then
                  call dualshape6tritilde
     &                 (xis,ets,xl2s,xsj2,xs2,shp2s,ns,
     &                 pslavdual,iflag)
                else
                  call dualshape6tritilde_lin
     &                 (xis,ets,xl2s,xsj2,xs2,shp2s,ns,
     &                 pslavdual,iflag)
                endif
              else
                call dualshape3tri(xis,ets,xl2s,xsj2,xs2,shp2s,ns,
     &               pslavdual,iflag)
              endif
              xim=pslavsurf(1,indexf+m)
              etm=pslavsurf(2,indexf+m)
!
              if(nopes.eq.8) then
                if((iflagdualquad.eq.2).or.(iflagdualquad.eq.4))then
                  call shape8qtilde(xim,etm,xl2s,xsj2m,xs2m,
     &                 shp2m,iflag)
                else
                  call shape4q(xim,etm,xl2s,xsj2m,xs2m,shp2m,
     &                 iflag)
                  dx=dsqrt(xsj2m(1)**2+xsj2m(2)**2+xsj2m(3)**2)
                  call shape8qtilde_lin(xim,etm,xl2s,xsj2m,
     &                 xs2m,shp2m,iflag)
                  dx=dsqrt(xsj2m(1)**2+xsj2m(2)**2+xsj2m(3)**2)
                endif
              elseif(nopes.eq.4) then
                call shape4q(xim,etm,xl2s,xsj2m,xs2m,shp2m,iflag)
              elseif(nopes.eq.6) then
                if((iflagdualquad.eq.2).or.(iflagdualquad.eq.4))then
                  call shape6tritilde(xim,etm,xl2s,xsj2m,
     &                 xs2m,shp2m,iflag)
                else
                  call shape6tritilde_lin(xim,etm,xl2s,xsj2m,
     &                 xs2m,shp2m,iflag)
                endif
              else
                call shape3tri(xim,etm,xl2s,xsj2m,xs2m,shp2m,iflag)
              endif
              if(m.eq.1)then
                do j=1,8
                  do jj=1,8
                    contr(8*(j-1)+jj)=0.0
                  enddo
                enddo
              endif
              dx=dsqrt(xsj2m(1)**2+xsj2m(2)**2+xsj2m(3)**2)
              do j=1,nopes
                ifs=getlocno(j,jfaces,nope)
                nodesf=kon(ipkon(nelems)+ifs)
                locs=j
                do jj=1,nopes
                  ifm=getlocno(jj,jfaces,nope)
                  nodem=kon(ipkon(nelems)+ifm)
                  locm=jj
                  contribution=shp2s(4,locs)*shp2m(4,locm)
     &                 *weight
     &                 *dx
                  contr(8*(j-1)+jj)=contr(8*(j-1)+jj)
     &                 +contribution
                  iscontr(8*(j-1)+jj)=nodesf
                  imcontr(8*(j-1)+jj)=nodem
                  contribution=0.d0
                enddo
              enddo
            enddo
            if(checkbiorthogonality) then
              write(*,*) 'checkbiorth.: contri,iscontr,imcontr',l
              do j=1,nopes
                write(*,105)(contr(8*(j-1)+jj),jj=1,nopes)
              enddo
            endif
          enddo
        endif
!
        checknorm=.false.
        if(checknorm)then
          do l=itiefac(1,i),itiefac(2,i)
            ifaces=islavsurf(1,l)
            nelems=int(ifaces/10)
            jfaces=ifaces-nelems*10
            call getnumberofnodes(nelems,jfaces,lakon,nope,
     &           nopes,idummy)
            mint2d=islavsurf(2,l+1)-islavsurf(2,l)
            if(mint2d.eq.0) cycle
            indexf=islavsurf(2,l)
!
            do j=1,nope
              konl(j)=kon(ipkon(nelems)+j)
            enddo
            do m=1,nopes
              do j=1,3
                ifac=getlocno(m,jfaces,nope)
                xl2s(j,m)=co(j,konl(ifac))+
     &               vold(j,konl(ifac))
              enddo
            enddo
            if((iflagdualquad.eq.1).or.(iflagdualquad.eq.3))then
              if(nopes.eq.6)then
                nopes=3
              else if(nopes.eq.8)then
                nopes=4
              endif
            endif
            do m=1,mint2d
              xis=pslavsurf(1,indexf+m)
              ets=pslavsurf(2,indexf+m)
              weight=pslavsurf(3,indexf+m)
              ns=l
              iflag=2
              if(nopes.eq.8) then
                call dualshape8qtilde(xis,ets,xl2s,xsj2,xs2,
     &               shp2s,ns,pslavdual,iflag)
              elseif(nopes.eq.4) then
                call dualshape4q(xis,ets,xl2s,xsj2,xs2,shp2s,ns,
     &               pslavdual,iflag)
              elseif(nopes.eq.6) then
                call dualshape6tritilde(xis,ets,xl2s,xsj2,xs2,
     &               shp2s,ns,pslavdual,iflag)
              else
                call dualshape3tri(xis,ets,xl2s,xsj2,xs2,shp2s,ns,
     &               pslavdual,iflag)
              endif
              if(m.eq.1)then
                do j=1,nopes
                  contr(j)=0.0
                enddo
              endif
              dx= dsqrt(xsj2(1)**2+xsj2(2)**2+xsj2(3)**2)
!
              do j=1,nopes
                ifs=getlocno(j,jfaces,nope)
                nodesf=kon(ipkon(nelems)+ifs)
                locs=j
                contribution=shp2s(4,locs)*
     &               weight
     &               *dx
                contr(j)=contr(j)
     &               +contribution
                iscontr(j)=nodesf
                contribution=0.d0
              enddo
            enddo
            write(*,*) 'checknorm: contri,iscontr,imcontr',l
            do j=1,nopes
              write(*,*)contr(j),iscontr(j)
            enddo
!
          enddo
        endif
      enddo
c      write(*,*) 'gendualcoeffs.f'
c      do j=1,64
c        write(*,*) j,pslavdual(j,1)
c      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
!****************************************************
 module string_routines
!****************************************************
 implicit none
 private
 public :: strip_comments
 contains
!****************************************************

	 function strip_comments(str,c) result(str2)
	 implicit none
	 character(len=*),intent(in) :: str
	 character(len=1),intent(in) :: c !comment character
	 character(len=len(str)) :: str2

	 integer :: i

	 i = index(str,c)
	 if (i>0) then
		str2 = str(1:i-1)
	 else
		str2 = str
	 end if

	 end function strip_comments

!****************************************************
 end module string_routines
!****************************************************

!****************************************************
 program main
!****************************************************
! Example use of strip_comments function
!****************************************************
 use string_routines, only: strip_comments
 implicit none

 write(*,*) strip_comments('apples, pears # and bananas', '#')
 write(*,*) strip_comments('apples, pears ; and bananas', ';')

!****************************************************
 end program main
!****************************************************

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
!##############################################################################
!# ****************************************************************************
!# <name> assemblytemplates </name>
!# ****************************************************************************
!#
!# <purpose>
!# Contains general template structures that are used to collect information
!# necessary for the assembly of matrices and vectors.
!#
!# The following routines can be found here:
!#
!# 1.) astmpl_createLevelInfoHier
!#     -> Creates a template assembly hierarchy
!#
!# 2.) astmpl_releaseLevelInfoHier
!#     -> Releases a template assembly hierarchy
!#
!# </purpose>
!##############################################################################

module assemblytemplates

  use fsystem
  use storage
  use boundary
  use triangulation
  use linearsystemscalar

  use spatialdiscretisation
  use timediscretisation
  use timescalehierarchy

  use meshhierarchy
  use fespacehierarchybase
  use fespacehierarchy
  use spacetimehierarchy

  implicit none

  private

!<types>

!<typeblock>

  ! A type block specifying all "static" information in space which are depending
  ! on a discretisation and a triangulation. Such static information can be
  ! precalculated and is valid until the mesh or the FE spaces change.
  type t_staticSpaceAsmTemplates

    ! An object for saving the triangulation on the domain
    type(t_triangulation), pointer :: p_rtriangulation => null()

    ! A scalar discretisation structure for the velocity space.
    type(t_spatialDiscretisation), pointer :: p_rdiscr => null()

    ! A scalar discretisation structure for the pressure space.
    type(t_spatialDiscretisation), pointer :: p_rdiscrPressure => null()

    ! A scalar discretisation structure that specifies how to generate
    ! the mass matrix in the velocity FEM space.
    ! May use a different cubature rule that the velocity discretisation
    type(t_spatialDiscretisation), pointer :: p_rdiscrMass => null()

    ! A scalar discretisation structure that specifies how to generate
    ! the mass matrix in the pressure FEM space.
    type(t_spatialDiscretisation), pointer :: p_rdiscrMassPressure => null()

    ! Cubature information structure for Mass-type matrices in the velocity space.
    type(t_scalarCubatureInfo) :: rcubatureInfoMass

    ! Cubature information structure for Laplace-type matrices.
    type(t_scalarCubatureInfo) :: rcubatureInfo

    ! Cubature information structure for pressure/divergence matrices.
    type(t_scalarCubatureInfo) :: rcubatureInfoDiv

    ! Cubature information structure for Mass-type matrices in the velocity space.
    type(t_scalarCubatureInfo) :: rcubatureInfoMassPressure

    ! Cubature information structure for RHS vectors, momentum equation
    type(t_scalarCubatureInfo) :: rcubatureInfoRHS

    ! Cubature information structure for RHS vectors, continuity equation
    type(t_scalarCubatureInfo) :: rcubatureInfoRHScontinuity

    ! A template FEM matrix that defines the structure of Laplace/Stokes/...
    ! matrices. The matrix contains only a stucture, no content.
    type(t_matrixScalar) :: rmatrixTemplateFEM

    ! A template FEM matrix that defines the structure of the offdiagonal
    ! matrices. Needed e.g. by Newton and may have different structure
    ! than the diagonal velocity matrices in case EOJ is activated.
    ! The matrix contains only a stucture, no content.
    type(t_matrixScalar) :: rmatrixTemplateFEMOffdiag

    ! A template FEM matrix that defines the structure of the pressure
    ! matrices. The matrix contains only a stucture, no content.
    type(t_matrixScalar) :: rmatrixTemplateFEMPressure

    ! A template FEM matrix that defines the structure of gradient
    ! matrices (B1/B2). The matrix contains only a stucture, no content.
    type(t_matrixScalar) :: rmatrixTemplateGradient

    ! A template FEM matrix that defines the structure of divergence
    ! matrices (D1/D2). The matrix contains only a stucture, no content.
    type(t_matrixScalar) :: rmatrixTemplateDivergence

    ! Precalculated Laplace matrix for that specific level
    type(t_matrixScalar) :: rmatrixLaplace

    ! Precalculated B1-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixB1

    ! Precalculated B2-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixB2

    ! Precalculated D1-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixD1

    ! Precalculated D2-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixD2

    ! Precalculated D1^T-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixD1T

    ! Precalculated D2-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixD2T

    ! Precalculated mass matrix
    type(t_matrixScalar) :: rmatrixMassPrimal

    ! Precalculated lumped mass matrix.
    type(t_matrixScalar) :: rmatrixMassLumpInt

    ! Precalculated mass matrix for the pressure space.
    type(t_matrixScalar) :: rmatrixMassPressure

    ! Precalculated mass; control space
    type(t_matrixScalar) :: rmatrixMassDistC

    ! Precalculated lumped mass; control space
    type(t_matrixScalar) :: rmatrixMassDistCLumped

    ! Precalculated mass matrix for the pressure space.
    ! Extended structure applicable for UMFPACK on pure Dirichlet problems.
    type(t_matrixScalar) :: rmatrixMassPressureExtStruc

    ! Precalculated lumped mass matrix for the pressure space.
    type(t_matrixScalar) :: rmatrixMassPressureLumpInt

    ! Precalculated EOJ matrix in the primal velocity space.
    ! Calculated with nu=1.
    type(t_matrixScalar) :: rmatrixEOJPrimal

    ! Precalculated EOJ matrix in the dual velocity space.
    ! Calculated with nu=1. May share the content with rmatrixEOJPrimal.
    type(t_matrixScalar) :: rmatrixEOJDual
  end type

!</typeblock>

  public :: t_staticSpaceAsmTemplates

!<typeblock>

  ! A hierarchy of t_staticSpaceAsmTemplates structures.
  type t_staticSpaceAsmHierarchy

    ! Number of levels in the hierarchy.
    integer :: nlevels = 0

    ! The level info structures on all levels.
    type(t_staticSpaceAsmTemplates), dimension(:), pointer :: p_RasmTemplList => null()

  end type

!</typeblock>

  public :: t_staticSpaceAsmHierarchy

!</types>

  ! Allocates memory for a level hierarchy consisting of nlevels levels
  public :: astmpl_createSpaceAsmHier

  ! Releases a level info hierarchy.
  public :: astmpl_releaseSpaceAsmHier

contains

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

!<subroutine>

  subroutine astmpl_createSpaceAsmHier (rhierarchy,nlevels)

!<description>
  ! Allocates memory for a level hierarchy consisting of nlevels levels
!</description>

!<input>
  ! Number of levels.
  integer :: nlevels
!</input>

!<output>
  ! A t_staticLevelInfoHierarchy to initialise.
  type(t_staticSpaceAsmHierarchy), intent(out) :: rhierarchy
!</output>

!</subroutine>

    ! Allocate memory.
    rhierarchy%nlevels = nlevels
    allocate(rhierarchy%p_RasmTemplList(nlevels))

  end subroutine

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

!<subroutine>

  subroutine astmpl_releaseSpaceAsmHier (rhierarchy)

!<description>
  ! Releases a level info hierarchy.
  !
  ! WARNING: Attached matrices are NOT automatically released!
!</description>

!<inputoutput>
  ! A t_staticLevelInfoHierarchy to clean up.
  type(t_staticSpaceAsmHierarchy), intent(inout) :: rhierarchy
!</inputoutput>

!</subroutine>

    ! Release memory.
    deallocate(rhierarchy%p_RasmTemplList)
    rhierarchy%nlevels = 0

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
module gauss_seidel

  use shared_data
  use boundary_conditions
  use diagnostics

  implicit none

  private

  public :: solve_const_Helmholtz, &
            solve_variable_elliptic

  contains

  subroutine solve_const_Helmholtz(phigs,f,alpha,beta,use_old_phi,tol)

    ! Solves the constand co-efficient Helmholtz equation
    ! (alpha - beta del**2) phi = f
    !
    ! For somple Poisson set alpha = 0 and beta = -1
    !
    ! Variable coefficients are handled separately

    real(num), dimension(-1:nx+2,-1:ny+2), intent(inout) :: phigs
    real(num), dimension(1:nx,1:ny), intent(in) :: f
    real(num), intent(in) :: alpha, beta
    logical, intent(in) :: use_old_phi
    real(num), intent(in) :: tol

    logical :: gsrb=.true. !hidden option

    real(num) :: L2, L2_old !norms
    real(num) :: L_phigs !laplacian of phigs

    integer :: maxir = -1 ! <0 for no max no of iterations
    integer :: ir = 0
    logical :: verbose=.false.

    print *, '*** begining GS relaxation solve.'
    print *, '*** this can take a while, use VERBOSE if you want to monitor stepping'

    if (.not. use_old_phi) then
      phi = 0.0_num
    else
      call phigs_bcs(phigs) !phi = 0.0_num
    endif

    L2_old = 1e6_num

    ir = 0
    do
      ir = ir + 1

      ! if not using redblack order
      if (.not. gsrb) then
        do iy = 1, ny
        do ix = 1, nx
          phigs(ix,iy) = phigs(ix+1,iy) + phigs(ix-1,iy) + phigs(ix,iy+1) + phigs(ix,iy-1)
          phigs(ix,iy) = (f(ix,iy) + beta * phigs(ix,iy) / dx**2) / &
            & (alpha + 4.0_num * beta / dx**2)

! hardcoded poisson equation
!          phi(ix,iy) = 0.25_num * ( &
!            & phi(ix+1,iy) + phi(ix-1,iy) + phi(ix,iy+1) + phi(ix,iy-1) &
!            - dx**2 * divu(ix,iy) )
        end do
        end do
      else !use red black
        ! odd iteration
        do iy = 1, ny
        do ix = 1, nx
          if (modulo(ix+iy,2) == 1) then
            phigs(ix,iy) = phigs(ix+1,iy) + phigs(ix-1,iy) + phigs(ix,iy+1) + phigs(ix,iy-1)
            phigs(ix,iy) = (f(ix,iy) + beta * phigs(ix,iy) / dx**2) / &
              & (alpha + 4.0_num * beta / dx**2)

!            phi(ix,iy) = 0.25_num * ( &
!              & phi(ix+1,iy) + phi(ix-1,iy) + phi(ix,iy+1) + phi(ix,iy-1) &
!              - dx**2 * divu(ix,iy) )
          endif
        end do
        end do
        ! even iteration
        do iy = 1, ny
        do ix = 1, nx
          if (modulo(ix+iy,2) == 0) then
            phigs(ix,iy) = phigs(ix+1,iy) + phigs(ix-1,iy) + phigs(ix,iy+1) + phigs(ix,iy-1)
            phigs(ix,iy) = (f(ix,iy) + beta * phigs(ix,iy) / dx**2) / &
              & (alpha + 4.0_num * beta / dx**2)
!            phi(ix,iy) = 0.25_num * ( &
!              & phi(ix+1,iy) + phi(ix-1,iy) + phi(ix,iy+1) + phi(ix,iy-1) &
!              - dx**2 * divu(ix,iy) )
          endif
        end do
        end do

      endif

      ! Apply periodic boundary conditions on phi's ghost cells

      call phigs_bcs(phigs)

      L2 = 0.0_num
      do iy = 1, ny
      do ix = 1, nx
        L_phigs = (phigs(ix+1,iy) - 2.0_num*phigs(ix,iy) + phigs(ix-1,iy)) &
          & / dx**2 + &
          & (phigs(ix,iy+1) - 2.0_num*phigs(ix,iy) + phigs(ix,iy-1)) / dy**2

        !L2 = L2 + abs(divu(ix,iy)-L_phi)**2
        L2 = L2 + abs(f(ix,iy)-(alpha*phigs(ix,iy)-beta*L_phigs))**2
      end do
      end do
      L2 = sqrt( L2 / REAL(nx*ny,num))

      if (verbose) &
        & print *, 'GS-iteration',ir,'complete. L2 is',L2,'|L2-L2_old| is',abs(L2-L2_old),'tol is',tol

     !exit conditions
     !if ((step >= nsteps .and. nsteps >= 0) .or. (L2 <= tol)) exit
     ! alt, exit if the difference between L2 and L2 prev is small - might
     ! indicate convergence
      if ((ir >= maxir .and. maxir >= 0) .or. (abs(L2-L2_old) <= tol)) exit
      L2_old = L2
    end do

    print *, '*** Gauss Seidel relaxation completed in',ir,'steps'
    print *, '*** L2 norm on numerical residuals',L2

  end subroutine solve_const_Helmholtz

  subroutine solve_variable_elliptic(phigs,f,eta,use_old_phi,tol)

    real(num), dimension(-1:nx+2,-1:ny+2), intent(inout) :: phigs
    real(num), dimension(1:nx,1:ny), intent(in) :: f
    real(num), dimension(0:nx+1,0:ny+1), intent(in) :: eta
            ! ^ eta should typically be rho plus its up-to-date +-1 ghost cells
    logical, intent(in) :: use_old_phi
    real(num), intent(in) :: tol

    logical :: gsrb=.true. !hidden option

    real(num) :: L2, L2_old !norms
    real(num) :: L_phigs !laplacian of phigs

    real(num) :: eta_ip, eta_im, eta_jp, eta_jm

    integer :: maxir = -1 ! <0 for no max no of iterations
    integer :: ir = 0
    logical :: verbose=.false.

    print *, '*** begining GS relaxation solve.'
    print *, '*** this can take a while, use VERBOSE if you want to monitor stepping'

    if (.not. use_old_phi) then
      phi = 0.0_num
    else
      call phigs_bcs(phigs) !phi = 0.0_num
    endif

    L2_old = 1e6_num

    ir = 0
    do
      ir = ir + 1

      ! if not using redblack order
      if (.not. gsrb) then
        do iy = 1, ny
        do ix = 1, nx
          !eta_ip actually also has the deltas built in
          eta_ip = 0.5_num * (eta(ix,iy)+eta(ix+1,iy)) / dx**2
          eta_im = 0.5_num * (eta(ix,iy)+eta(ix-1,iy)) / dx**2
          eta_jp = 0.5_num * (eta(ix,iy)+eta(ix,iy+1)) / dy**2
          eta_jm = 0.5_num * (eta(ix,iy)+eta(ix,iy-1)) / dy**2

          phigs(ix,iy) = eta_ip * phigs(ix+1,iy) + eta_im*phigs(ix-1,iy) &
            & + eta_jp * phigs(ix,iy+1) + eta_jm * phigs(ix,iy-1) - f(ix,iy)
          phigs(ix,iy) = phigs(ix,iy) / ( eta_ip + eta_im + eta_jp + eta_jm )
        end do
        end do
      else !use red black
        ! odd iteration
        do iy = 1, ny
        do ix = 1, nx
          if (modulo(ix+iy,2) == 1) then
            !eta_ip actually also has the deltas built in
            eta_ip = 0.5_num * (eta(ix,iy)+eta(ix+1,iy)) / dx**2
            eta_im = 0.5_num * (eta(ix,iy)+eta(ix-1,iy)) / dx**2
            eta_jp = 0.5_num * (eta(ix,iy)+eta(ix,iy+1)) / dy**2
            eta_jm = 0.5_num * (eta(ix,iy)+eta(ix,iy-1)) / dy**2

            phigs(ix,iy) = eta_ip * phigs(ix+1,iy) + eta_im*phigs(ix-1,iy) &
              & + eta_jp * phigs(ix,iy+1) + eta_jm * phigs(ix,iy-1) - f(ix,iy)
            phigs(ix,iy) = phigs(ix,iy) / ( eta_ip + eta_im + eta_jp + eta_jm )
          endif
        end do
        end do
        ! even iteration
        do iy = 1, ny
        do ix = 1, nx
          if (modulo(ix+iy,2) == 0) then
            !eta_ip actually also has the deltas built in
            eta_ip = 0.5_num * (eta(ix,iy)+eta(ix+1,iy)) / dx**2
            eta_im = 0.5_num * (eta(ix,iy)+eta(ix-1,iy)) / dx**2
            eta_jp = 0.5_num * (eta(ix,iy)+eta(ix,iy+1)) / dy**2
            eta_jm = 0.5_num * (eta(ix,iy)+eta(ix,iy-1)) / dy**2

            phigs(ix,iy) = eta_ip * phigs(ix+1,iy) + eta_im*phigs(ix-1,iy) &
              & + eta_jp * phigs(ix,iy+1) + eta_jm * phigs(ix,iy-1) - f(ix,iy)
            phigs(ix,iy) = phigs(ix,iy) / ( eta_ip + eta_im + eta_jp + eta_jm )
          endif
        end do
        end do

      endif

      ! Apply periodic boundary conditions on phi's ghost cells

      call phigs_bcs(phigs)

      L2 = 0.0_num
      do iy = 1, ny
      do ix = 1, nx

        eta_ip = 0.5_num * (eta(ix,iy)+eta(ix+1,iy))
        eta_im = 0.5_num * (eta(ix,iy)+eta(ix-1,iy))
        eta_jp = 0.5_num * (eta(ix,iy)+eta(ix,iy+1))
        eta_jm = 0.5_num * (eta(ix,iy)+eta(ix,iy-1))

        L_phigs = ( eta_ip * (phigs(ix+1,iy)-phigs(ix,iy)) &
              &    -  eta_im * (phigs(ix,iy)-phigs(ix-1,iy)) ) / dx**2 + &
              & ( eta_jp * (phigs(ix,iy+1)-phigs(ix,iy)) &
              &    -  eta_jm * (phigs(ix,iy)-phigs(ix,iy-1)) ) / dy**2

        !residual
        L2 = L2 + abs(f(ix,iy)-L_phigs)**2
      end do
      end do
      L2 = sqrt( L2 / REAL(nx*ny,num))

      if (verbose) &
        & print *, 'GS-iteration',ir,'complete. L2 is',L2,'|L2-L2_old| is',abs(L2-L2_old),'tol is',tol

     !exit conditions
     !if ((step >= nsteps .and. nsteps >= 0) .or. (L2 <= tol)) exit
     ! alt, exit if the difference between L2 and L2 prev is small - might
     ! indicate convergence
      if ((ir >= maxir .and. maxir >= 0) .or. (abs(L2-L2_old) <= tol)) exit
      L2_old = L2
    end do

    print *, '*** Gauss Seidel relaxation completed in',ir,'steps'
    print *, '*** L2 norm on numerical residuals',L2

  end subroutine solve_variable_elliptic

  subroutine phigs_bcs(phigs)

    real(num), dimension(-1:nx+2,-1:ny+2), intent(inout) :: phigs

    if (bc_xmin == periodic) then
      phigs(0,:) = phigs(nx,:)
      phigs(-1,:) = phigs(nx-1,:)
    endif
    if (bc_xmax == periodic) then
      phigs(nx+1,:) = phigs(1,:)
      phigs(nx+2,:) = phigs(2,:)
    endif
    if (bc_ymin == periodic) then
      phigs(:,0) = phigs(:,ny)
      phigs(:,-1) = phigs(:,ny-1)
    endif
    if (bc_ymax == periodic) then
      phigs(:,ny+1) = phigs(:,1)
      phigs(:,ny+2) = phigs(:,2)
    endif

    ! Zero gradient

    if (bc_xmin == zero_gradient) then
      phigs(0,:) = phigs(1,:)
      phigs(-1,:) = phigs(2,:)
    endif
    if (bc_xmax == zero_gradient) then
      phigs(nx+1,:) = phigs(nx,:)
      phigs(nx+2,:) = phigs(nx-1,:)
    endif
    if (bc_ymin == zero_gradient) then
      phigs(:,0) = phigs(:,1)
      phigs(:,-1) = phigs(:,2)
    endif
    if (bc_ymax == zero_gradient) then
      phigs(:,ny+1) = phigs(:,ny)
      phigs(:,ny+2) = phigs(:,ny-1)
    endif

    ! No slip

!!    if (bc_xmin == no_slip) then
!!      phigs(0,:) = -phigs(1,:)
!!      phigs(-1,:) = -phigs(2,:)
!!    endif
!!    if (bc_xmax == no_slip) then
!!      phigs(nx+1,:) = -phigs(nx,:)
!!      phigs(nx+2,:) = -phigs(nx-1,:)
!!    endif
!!    if (bc_ymin == no_slip) then
!!      phigs(:,0) = -phigs(:,1)
!!      phigs(:,-1) = -phigs(:,2)
!!    endif
!!    if (bc_ymax == no_slip) then
!!      phigs(:,ny+1) = -phigs(:,ny)
!!      phigs(:,ny+2) = -phigs(:,ny-1)
!!    endif
    if (bc_xmin == no_slip) then
      phigs(0,:) = phigs(1,:)
      phigs(-1,:) = phigs(2,:)
    endif
    if (bc_xmax == no_slip) then
      phigs(nx+1,:) = phigs(nx,:)
      phigs(nx+2,:) = phigs(nx-1,:)
    endif
    if (bc_ymin == no_slip) then
      phigs(:,0) = phigs(:,1)
      phigs(:,-1) = phigs(:,2)
    endif
    if (bc_ymax == no_slip) then
      phigs(:,ny+1) = phigs(:,ny)
      phigs(:,ny+2) = phigs(:,ny-1)
    endif

!    !overwrite with experimental no_slip
    ! seems to give large (>0 !) for the grav problem

!    if (bc_ymin == no_slip) then
!      phigs(:,0) =  phigs(:,1) - (phigs(:,2)-phigs(:,1))
!      phigs(:,-1) = phigs(:,1) - (phigs(:,2)-phigs(:,1))*2.0_num
!    endif
!    if (bc_ymax == no_slip) then
!      phigs(:,ny+1) = phigs(:,ny) + (phigs(:,ny) - phigs(:,ny-1))
!      phigs(:,ny+2) = phigs(:,ny) + (phigs(:,ny) - phigs(:,ny-1))*2.0_num
!    endif

    !duplicate of no_slip for driven
    if (bc_ymax == driven) then
      phigs(:,ny+1) = phigs(:,ny)
      phigs(:,ny+2) = phigs(:,ny-1)
    endif

  end subroutine phigs_bcs

end module gauss_seidel

max_stars_repo_name: JOThurgood/SimpleCFD

max_issues_repo_name: JOThurgood/SimpleCFD

max_forks_repo_name: JOThurgood/SimpleCFD


text:
! { dg-do run }
! Fix for PR29699 - see below for details.
!
! Contributed by Tobias Burnus  <burnus@gcc.gnu.org>
!
PROGRAM vocabulary_word_count

  IMPLICIT NONE
  TYPE VARYING_STRING
    CHARACTER,DIMENSION(:),ALLOCATABLE :: chars
  ENDTYPE VARYING_STRING

  INTEGER :: list_size=200

  call extend_lists2

CONTAINS

! First the original problem: vocab_swap not being referenced caused
! an ICE because default initialization is used, which results in a
! call to gfc_conv_variable, which calls gfc_get_symbol_decl.

  SUBROUTINE extend_lists1
    type(VARYING_STRING),DIMENSION(list_size) :: vocab_swap
  ENDSUBROUTINE extend_lists1

! Curing this then uncovered two more problems: If vocab_swap were
! actually referenced, an ICE occurred in the gimplifier because
! the declaration for this automatic array is presented as a
! pointer to the array, rather than the array. Curing this allows
! the code to compile but it bombed out at run time because the
! malloc/free occurred in the wrong order with respect to the
! nullify/deallocate of the allocatable components.

  SUBROUTINE extend_lists2
    type(VARYING_STRING),DIMENSION(list_size) :: vocab_swap
    allocate (vocab_swap(1)%chars(10))
    if (.not.allocated(vocab_swap(1)%chars)) call abort ()
    if (allocated(vocab_swap(10)%chars)) call abort ()
  ENDSUBROUTINE extend_lists2

ENDPROGRAM vocabulary_word_count

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
*
*     -------------------------------------------------------------
*      N O N R E L A T 4 1
*     -------------------------------------------------------------
*                                                                  *
*     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
*     OF TWO PARTICLE OPERATOR IN CASE :       N'1 = N1 + 1        *
*                                              N'2 = N2 + 1        *
*                                              N'3 = N3 - 2,       *
*     WHEN IREZ = 1   ........................
*                                              N'1 = N1 - 1        *
*                                              N'2 = N2 - 1        *
*                                              N'3 = N3 + 2,       *
*     WHEN IREZ = 2   ........................
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Vanderbilt University,  Nashville            February 1994   *
*
      SUBROUTINE NONRELAT41(IA,IB,IC,IREZ,IIA,IIB,IIC,IID)
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
      LOGICAL RECOUPLS0,CALCULATION
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/MEDEFN/IHSH,NJ(16),LJ(16),NOSH1(16),NOSH2(16),J1QN1(31,3),
     :     J1QN2(31,3),IJFUL(16)
      COMMON/DIAGNL/IDIAG,JA,JB
      COMMON /OPERAT/ ICOLOM,ISOTOP,IORBORB
      COMMON/TRK/BD1(3),BD2(3),BK1(3),BK2(3),
     *ID1(7),ID2(7),IK1(7),IK2(7)
      COMMON/TRK2/BD3(3),BD4(3),BK3(3),BK4(3),
     *ID3(7),ID4(7),IK3(7),IK4(7)
      DIMENSION PMS(2),IATS(2)
      IF(IHSH.LE.2)RETURN
      CALL EILE(IA,IB,IC,IAA,IBB,ICC)
      IF(.NOT.RECOUPLS0(1,IAA,ICC,IBB,IBB,0))RETURN
      IF(.NOT.RECOUPLS0(2,IAA,ICC,IBB,IBB,2))RETURN
      IF(.NOT.RECOUPLS0(3,IAA,ICC,IBB,IBB,2))RETURN
      IF(IREZ.EQ.1) THEN
        QM1=-HALF
        QM2=HALF
      ELSE
        QM1=HALF
        QM2=-HALF
      ENDIF
      A1=ZERO
      A2=ZERO
      CALL HIBFF(IA,IB,IC,IA,3)
C
C     CASES 3312   + + - -        TRANSFORM TO  1233   - - + +
C           3321                                1233
C                                                    (IREZ = 1)
C     OR
C     CASES 1233   + + - -        TRANSFORM TO  1233   + + - -
C           2133                                1233
C                                                    (IREZ = 2)
      PMS(1)=ZERO
      PMS(2)=ZERO
      DO 1 I13=1,2
        JS12=I13-1
        CALL RECOUPLS3(3,IA,IB,IC,1,1,JS12*2,0,IAT,REC)
        IATS(I13)=IAT
        IF(IATS(I13).NE.0) THEN
          CALL RECOUPLS3(3,IA,IB,IC,1,1,JS12*2,1,IAT,RECS)
          PMS(I13)=RECS
        ENDIF
    1 CONTINUE
      IF((IATS(1)+IATS(2)).EQ.0)RETURN
      LA=IJFUL(IIA)
      LB=IJFUL(IIB)
      LC=IJFUL(IIC)
      LD=IJFUL(IID)
      LIA2=LJ(IA)*2
      LIB2=LJ(IB)*2
      LIC2=LJ(IC)*2
      IP1=ITREXG(LJ(IB),LJ(IA),LJ(IC),LJ(IC),IKK)+1
      IF(IKK.LE.0)RETURN
      IG1=IP1+IKK-1
      IP2=ITREXG(LJ(IC),LJ(IA),LJ(IB),LJ(IC),IKK)+1
      IF(IKK.LE.0)RETURN
      IG2=IP2+IKK-1
      IFAZP=JFAZE(IC,IA,IB,IC)
      IF(IA.LT.IB) THEN
         IAA=IA
         IBB=IB
      ELSE
         IAA=IB
         IBB=IA
      ENDIF
      NN=0
      IB1=IBB-1
      DO 5 II=IAA,IB1
        NN=NOSH1(II)+NN
    5 CONTINUE
      IF((NN/2)*2.EQ.NN)IFAZP=-IFAZP
      DO 2 I2=IP2,IG2
        KL=I2-1
	IF(CALCULATION(KL)) THEN
          IF((ICOLOM+ISOTOP).EQ.1)
     :           CALL COULOMBLS(LJ(IA),LJ(IB),LJ(IC),LJ(IC),KL,A1)
          IF(IORBORB.EQ.1)
     :           CALL ORBITORBIT(LJ(IA),LJ(IB),LJ(IC),LJ(IC),KL+1,A2)
          IF((DABS(A1)+DABS(A2)).GT.EPS) THEN
            AB=ZERO
            DO 3 I3=IP1,IG1
              L12=I3-1
              CALL RECOUPLS3(2,IA,IB,IC,LIA2,LIB2,L12*2,0,IAT,REC)
              IF(IAT.NE.0) THEN
                IF(IXJTIK(LIC2,LIB2,KL*2,LIA2,LIC2,L12*2).NE.0) THEN
                  AC=ZERO
                  DO 4 I13=1,2
		    IF(IATS(I13).NE.0) THEN
                      JS12=I13-1
		      IFAZ=L12+JS12
	              IF((IFAZ/2)*2.EQ.IFAZ) THEN
                        CALL A1A2W3LS(IK1,IK2,IK3,BK1,BK2,BK3,
     :          ID1,ID2,ID3,BD1,BD2,BD3,L12,JS12,QM1,QM1,QM2,QM2,AA)
                        AA=AA*PMS(I13)*DSQRT(DBLE(2*JS12+1))
		        AC=AC+AA
	              ENDIF
	            ENDIF
    4             CONTINUE
                  CALL RECOUPLS3(2,IA,IB,IC,LIA2,LIB2,L12*2,1,IAT,RECL)
                  CALL SIXJ(LIC2,LIB2,KL*2,LIA2,LIC2,L12*2,0,SI)
                  AA=AC*RECL*SI*DSQRT(DBLE(2*L12+1))
                  IFAZ=IK1(3)+IK3(3)+KL+L12
	          IF(IREZ.EQ.2) IFAZ=IK2(3)+IK3(3)+KL+L12
                  IF((IFAZ/2)*2.NE.IFAZ)AA=-AA
                  AB=AB+AA
	        ENDIF
	      ENDIF
    3       CONTINUE
            IF(DABS(AB).GT.EPS) THEN
              AB=-AB*DBLE(IFAZP)
	      ABB=AB
              IF((ICOLOM+ISOTOP).EQ.1) THEN
                IF(DABS(A1).GT.EPS) THEN
                  AB=A1*AB
                  CALL SAVENON(3,AB,KL,LA,LB,LC,LD,JA,JB,0)
                ENDIF
              ENDIF
C Orbit 3312
              IF(IORBORB.EQ.1) THEN
                IF(DABS(A2).GT.EPS) THEN
	          KLL=KL-1
                  ABB=A2*ABB*HALF/DSQRT(DBLE(2*KL+1))
                  CALL SAVENON(9,ABB,KLL,LA,LB,LC,LD,JA,JB,0)
                  CALL SAVENON(9,ABB,KLL,LB,LA,LD,LC,JA,JB,0)
C                 WRITE(79,656) ABB,KLL,LJ(IA),LJ(IB),JA,JB
  656 FORMAT(1X,'1221 O-O','AA=',F17.7,'K=',I3,'LA=',2I3,'JA=',I4,
     :'JB=',I4)
                END IF
              END IF
            ENDIF
          ENDIF
	ENDIF
    2 CONTINUE
      RETURN
      END


text:
module m
  integer::i=3
contains
  subroutine s
    print '(a)','hello'
  end subroutine s
end module m

program t
  use m,only:s
  integer::i=77
  print '(i0)',i
  call s
end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
      implicit double precision (a-h,o-z)
c
c             make SSY boundary displacements for K and T loading.
c             also scale existing coords to a different radius
c             if desired.
c
c
      dimension coord(3,5000), node_list(5000)
      double precision new_radius
      integer elem, outfile
      character * 80 fname, text * 1
c
      write(*,*) ' '
      write(*,*) '>> Compute applied displacements matching K & T'
      write(*,*) ' '
      write(*,*) ' > Coordinate file name: '
      read(*,fmt='(a80)') fname
      infile = 10
      open(unit=infile,file=fname,status='old')
      write(*,*)
     &   ' > Final K(I), T/sig-o, sig-o, number of load steps: '
      read(*,*) final_k, t_ratio, sig_o, num_steps
      applied_k = final_k / real(num_steps)
      applied_t = sig_o * t_ratio / real(num_steps)
      write(*,*) ' > Start radius, new scaled radius: '
      read(*,*) start_radius, new_radius
c
      read(infile,*) nnode
      do node = 1, nnode
       read(infile,*) ii, x, y, z
       x = x * (new_radius/start_radius)
       y = y * (new_radius/start_radius)
       coord(1,node) = x
       coord(2,node) = y
       coord(3,node) = z
      end do
      write(*,*) ' >> data file read'
c
c             loop to find all nodes on the outer most surface, i.e.,
c             those with radius.
c
      do node = 1, 5000
       node_list(node) = 0
      end do
c
      tolerance = 0.01
      node_count = 0
      do node = 1, nnode
       x = coord(1,node)
       y = coord(2,node)
       r = abs( sqrt( x*x + y*y ) - new_radius )
       if ( r .le. tolerance ) then
          node_count = node_count + 1
          node_list(node) = 1
       end if
      end do
      write(*,*) '>> Number of nodes on boundary: ',
     &          node_count
c
c             nodes to process have been marked in list.
c
      rnu = 0.3
      e = 30000.0
      pi = 3.14159
      open(unit=10,file='constraints',status='unknown')
      open(unit=11,file='new_coords',status='unknown')
      outfile = 10
      do node = 1, nnode
       x = coord(1,node)
       y = coord(2,node)
       if ( node_list(node) .eq. 1 ) then
            x = coord(1,node)
            y = coord(2,node)
            r = sqrt(x*x+y*y)
            theta = atan2( y, x )
            ct2 = cos(theta*0.5)
            st2 = sin(theta*0.5)
            ct  = cos(theta)
            st  = sin(theta)
            term1 = (3.0-4.0*rnu-ct) * ct2
            term2 = (3.0-4.0*rnu-ct) * st2
            term3 = sqrt( r * 0.5 / pi )
            term4 = (1.0 + rnu ) * applied_k / e
            u_k = term4 * term3 * term1
            v_k = term4 * term3 * term2
            u_t = (1.0-rnu*rnu) * applied_t * ct * r / e
            v_t = -rnu *(1.0+rnu) * applied_t * st * r / e
            u_total = u_k + u_t
            v_total = v_k + v_t
            write(outfile,1000) node, u_total, v_total
        end if
       write(11,1100) node, x, y, coord(3,node)
      end do
c
      close(unit=10)
      close(unit=11)
      write(*,*) '.. new coordinates is file: coords'
      write(*,*) '.. K and T constraints in file: constraints'
c
 1000 format(4x,i5,' u  ',e16.8,' v ',e16.8)
 1100 format(2x,i5,3e16.8)
      stop
      end