text:
      MODULE m_tetcon
      use m_juDFT
!
! This subroutine constructs the tetrahedra for the
! Brillouin zone integration
!
      CONTAINS
      SUBROUTINE tetcon(
     >                  iofile,ibfile,mkpt,ndiv3,
     >                  nkpt,omega,kvc,nsym,
     <                  nt,voltet,ntetra)

      IMPLICIT NONE
c
      INTEGER, INTENT (IN) :: iofile,ibfile
      INTEGER, INTENT (IN) :: mkpt,ndiv3,nkpt,nsym
      REAL,    INTENT (IN) :: omega
      REAL,    INTENT (IN) :: kvc(3,mkpt)
      INTEGER, INTENT (OUT) :: nt
      INTEGER, INTENT (OUT) :: ntetra(4,ndiv3)
      REAL,    INTENT (OUT) :: voltet(ndiv3)

      INTEGER i,it,j,ic,icom,ik,nkp,jj,nnt,nsid,nkq
      INTEGER i2,i3,i4,nk2,k,l,kk,is1,is2,is3,js1,js2,js3
      INTEGER nside(4),nside2(4),nkadd(mkpt)
      REAL dnorm,pi,vav,vsq,sav,ssq,dm,dist,dl,vect
      REAL vol,vt,length,minlen,eps,eps1
      REAL vmin,vmax,smin,smax,volnew,sum,sss
      REAL ddist(mkpt),kcorn(3,4),norm(3),d(3,16),dnm(3),cn(3)
C
C----->  Intrinsic Functions
C
       INTRINSIC abs,max,min,sqrt
C
C*************************************************************************
C
C     All commons have been removed from the historical program and re-
C     placed by direct subroutine calls.  Those statements immediately
C     below preceded by stars your removed prior to my changes
C
C                                          June 9, 1992
C                                                   Fred Hutson
C
C*************************************************************************
C
      data eps/1e-10/
      data eps1/1e-5/
      pi = 4.0*atan(1.0)
C
C CONSTRUCT THE FIRST TETRAHEDRON
C
      ntetra(1,1)=1
C
C FIND K-POINT NEAREST TO 1
C
      dm=10.0**9
      i2=0
      do 100 i=2,nkpt
      dist=0.0
      do 110 icom=1,3
      dist=dist+(kvc(icom,1)-kvc(icom,i))**2
  110 continue
      if ( dist.gt.dm ) goto 100
      i2=i
      dm=dist
  100 continue
      IF ( i2==0 )  CALL juDFT_error(" tetcon1 ",calledby ="tetcon")
      dnorm=sqrt(dm)
      ntetra(2,1)=i2
C
C FIND POINT NEAREST TO (1+I2)/2 , NOT ON THE LINE
C CONNECTING 1 AND I2
C
      dm=10.0**9
      i3=0
      do 200 i=2,nkpt
      if ( i.eq.i2 ) goto 200
      dl=0.0
      dist=0.0
      do 210 icom=1,3
      vect=kvc(icom,i)-0.5*(kvc(icom,1)+kvc(icom,i2))
      dist=dist+vect*vect
      dl=dl+vect*(kvc(icom,1)-kvc(icom,i2))
  210 continue
      dl=dl/dnorm
      if ( abs(dist-dl*dl).lt.0.01*dist ) goto 200
      if ( dist.gt.dm ) goto 200
      dm=dist
      i3=i
  200 continue
      IF ( i3==0 )  CALL juDFT_error(" tetcon2 ",calledby ="tetcon")
      ntetra(3,1)=i3
C
C FIND POINT NEAREST TO (1+I2+I3)/3 , NOT IN THE
C PLANE (1,I2,I3)
C
      cn(1)=(kvc(2,1)-kvc(2,i2))*(kvc(3,i2)-kvc(3,i3))-
     $      (kvc(3,1)-kvc(3,i2))*(kvc(2,i2)-kvc(2,i3))
      cn(2)=(kvc(3,1)-kvc(3,i2))*(kvc(1,i2)-kvc(1,i3))-
     $      (kvc(1,1)-kvc(1,i2))*(kvc(3,i2)-kvc(3,i3))
      cn(3)=(kvc(1,1)-kvc(1,i2))*(kvc(2,i2)-kvc(2,i3))-
     $      (kvc(2,1)-kvc(2,i2))*(kvc(1,i2)-kvc(1,i3))
      dnorm=0.0
      do 300 icom=1,3
      dnorm=dnorm+cn(icom)**2
  300 continue
      dnorm=1.0/sqrt(dnorm)
      do 310 icom=1,3
      cn(icom)=cn(icom)*dnorm
  310 continue
      i4=0
      dm=10.0**9
      do 400 i=2,nkpt
      if ( (i.eq.i2).or.(i.eq.i3) ) goto 400
      dist=0.0
      dl=0.0
      do 410 icom=1,3
      vect=kvc(icom,i)-(kvc(icom,1)+kvc(icom,i2)+
     $     kvc(icom,i3))/3.0
      dist=dist+vect**2
      dl=dl+vect*cn(icom)
  410 continue
      if ( dl*dl.lt.0.01*dist ) goto 400
      if ( dist.gt.dm ) goto 400
      i4=i
      dm=dist
      vt=dl/(dnorm*6.0)
  400 continue
      IF ( i4==0 )  CALL juDFT_error(" tetcon3 ",calledby ="tetcon")
      ntetra(4,1)=i4
      voltet(1)=abs(vt)
C
C ENTER LOOP FOR CONSTRUCTION OF TETRAHEDRA:
C
      nt=1
      it=0
 1000 continue
      it=it+1
C
C CHECK THE SIDES OF TETRAHEDRON IT
C
      do 1100 j=1,4
C
C CHECK SIDE OPPOSITE TO CORNER J:
C
      ic=0
      do 1200 i=1,4
      if ( i.eq.j ) goto 1200
      ic=ic+1
      nside(ic)=ntetra(i,it)
      do 1300 icom=1,3
      kcorn(icom,ic)=kvc(icom,ntetra(i,it))
 1300 continue
 1200 continue
      nside(4)=ntetra(j,it)
      is1=min(nside(1),nside(2),nside(3))
      is3=max(nside(1),nside(2),nside(3))
      is2=nside(1)+nside(2)+nside(3)-is1-is3
C
C CHECK IF THERE IS ALREADY A TETRAHEDRON CONNECTED
C TO THIS SIDE:
C
      do nnt=1,nt
        do 1310 nsid=1,4
          if ( nnt.eq.it ) goto 1310
          ic=0
          do 1320 i=1,4
            if ( i.eq.nsid ) goto 1320
            ic=ic+1
            nside2(ic)=ntetra(i,nnt)
 1320     continue
          js1=min(nside2(1),nside2(2),nside2(3))
          js3=max(nside2(1),nside2(2),nside2(3))
          js2=nside2(1)+nside2(2)+nside2(3)-js1-js3
          if ( (is1.eq.js1).and.(is2.eq.js2).and.(is3.eq.js3) )
     $      goto 1100
 1310   continue
      end do
C
C CONSTRUCT THE OUTWARD NORMAL ON THIS SIDE
C
      norm(1)=(kcorn(2,2)-kcorn(2,1))*(kcorn(3,3)-kcorn(3,1))-
     $        (kcorn(3,2)-kcorn(3,1))*(kcorn(2,3)-kcorn(2,1))
      norm(2)=(kcorn(3,2)-kcorn(3,1))*(kcorn(1,3)-kcorn(1,1))-
     $        (kcorn(1,2)-kcorn(1,1))*(kcorn(3,3)-kcorn(3,1))
      norm(3)=(kcorn(1,2)-kcorn(1,1))*(kcorn(2,3)-kcorn(2,1))-
     $        (kcorn(2,2)-kcorn(2,1))*(kcorn(1,3)-kcorn(1,1))
      vol=0.0
      do 1400 i=1,3
      vol=vol+norm(i)*(kvc(i,ntetra(j,it))-kcorn(i,1))
 1400 continue
      vol=vol/6.0
      if ( vol.lt.0.0 ) goto 1500
      do 1600 i=1,3
      norm(i)=-norm(i)
 1600 continue
 1500 continue
      vol=abs(vol)
C
C STORE THE K-POINT ADDRESS IN NKADD ARRAY ACCORDING TO THE
C ORDER OF DISTANCE BETWEEN THE K-POINT AND THIS FACE.
C
      do 1800 nkp=1,nkpt
      length=0.0
      do 1850 i=1,3
      length=length+(kvc(i,nkp)
     &     -(kcorn(i,1)+kcorn(i,2)+kcorn(i,3))/3.0)**2
 1850 continue
      ddist(nkp)=length
      nkadd(nkp)=nkp
 1800 continue
      do 1900 nkp=1,nkpt-1
      minlen=ddist(nkp)
      ik=nkp
      do 1950 nkq=nkp+1,nkpt
      if(ddist(nkq).gt.minlen-eps) goto 1950
      minlen=ddist(nkq)
      ik=nkq
 1950 continue
      ddist(ik)=ddist(nkp)
      k=nkadd(nkp)
      nkadd(nkp)=nkadd(ik)
      nkadd(ik)=k
 1900 continue
C
C CONSTRUCT A TETRAHEDRON WHICH CONNECT THIS FACE TO
C A K-POINT.
C
      ik=0
      do 2000 nkp=1,nkpt
      do 2010 i=1,3
      if ( nkadd(nkp).eq.nside(i) ) goto 2000
 2010 continue
      do 2050 i=1,3
      kcorn(i,4)=kvc(i,nkadd(nkp))
 2050 continue
      nside(4)=nkadd(nkp)
      vt=0.0
      do 2100 i=1,3
      vt=vt+norm(i)*(kcorn(i,4)-kcorn(i,1))
 2100 continue
      vt=vt/6.0
C
C REJECT POINT NKP IF IT IS ON THE WRONG SIDE .
C
      if ( vt.lt.eps*vol ) goto 2000
C
C CHECK IF THIS TETRAHEDRON INTERSECTS AN EXISTING ONE
C
      do 3100 nk2=1,nt
      if(nk2.eq.it) goto 3100
C
C FIRST CHECK IF TETRAHEDRON NK2 IS ON THE DANGEROUS SIDE
C
      do 3150 i=1,4
      sum=0.0
      do 3160 jj=1,3
      sum=sum+norm(jj)*(kvc(jj,ntetra(i,nk2))-kcorn(jj,1))
 3160 continue
      if ( sum.gt.eps ) goto 3170
 3150 continue
      goto 3100
 3170 continue
C
C TETRAHEDRON NK2 IS POTENTIALLY SUSPECT
C
      l=0
      do i=1,4
        do 3200 jj=1,4
          if(ntetra(i,nk2).eq.nside(jj)) goto 3200
          sum=0.0
          do icom=1,3
            cn(icom)=kvc(icom,ntetra(i,nk2))-kcorn(icom,jj)
            sum=sum+cn(icom)*cn(icom)
          end do
          do icom=1,3
            cn(icom)=cn(icom)/sqrt(sum)
          end do
          do k=1,l
            sum=0.0
            do icom=1,3
              sum=sum+cn(icom)*d(icom,k)
            end do
            if((1.0-sum).lt.eps) goto 3200
C
C EXCLUDE THE VECTOR WHICH HAS THE SAME DIRECTION AS AN EXISTING ONE.
C
          end do
          l=l+1
          do icom=1,3
            d(icom,l)=cn(icom)
          end do
 3200   continue
      end do
      IF(l<4)  CALL juDFT_error(" tetcon9 ",calledby ="tetcon")
C
C HERE, WE HAVE A SET OF D-VECTORS WHICH CONNECT THE CORNER POINTS
C OF ONE TETRAHEDRON NK2 WITH THOSE OF THE CURRENT TETRAHEDRON TO
C BE CHECKED.
C FIND A PLANE OF WHICH ALL D-VECTORS EXIST IN ONE SIDE.
C
      do 3510 i=1,l-1
        do 3500 jj=i+1,l
          dnm(1)=d(2,i)*d(3,jj)-d(3,i)*d(2,jj)
          dnm(2)=d(3,i)*d(1,jj)-d(1,i)*d(3,jj)
          dnm(3)=d(1,i)*d(2,jj)-d(2,i)*d(1,jj)
C
C DNM IS THE NORMAL VECTOR TO THE PLANE GIVEN BY
C I-TH AND JJ-TH D-VECTORS.
C
          if((abs(dnm(1)).lt.eps).and.(abs(dnm(2)).lt.eps)
     &       .and.(abs(dnm(3)).lt.eps)) goto 3500
          do 3600 k=1,l
            if((k.eq.i).or.(k.eq.jj)) goto 3600
            sum=0.0
            do icom=1,3
              sum=sum+d(icom,k)*dnm(icom)
            end do
            if(abs(sum).lt.eps) goto 3600
            do 3700 kk=1,l
              if((kk.eq.i).or.(kk.eq.jj).or.(kk.eq.k)) goto 3700
              sss=0.0
              do icom=1,3
                sss=sss+d(icom,kk)*dnm(icom)
              end do
              if(abs(sss).lt.eps) goto 3700
              if(sss*sum.lt.0.0) goto 3500
C
C IF K-TH AND KK-TH D-VECTORS EXIST IN OPPOSITE SIDE WITH
C RESPECT TO THE PLANE GIVEN BY I-TH AND JJ-TH D-VECTORS,
C WE WILL ATTEMPT THE NEXT PLANE. (GOTO 3500)
C
 3700       continue
            goto 3100
C
C WE SUCCEED WHEN THE PLANE SATISFIES THE CONDITION,
C I.E., THE CURRENT TETRAHEDRON DOES NOT INTERSECTS NK2-TH ONE:
C WE TRY TO CHECK THE NEXT TETRAHEDRON. (GOTO 3100)
C
 3600     continue
 3500   continue
 3510 continue
      goto 2000
C
C HERE, A TETRAHEDRON MADE OF J-TH FACE AND NKADD(NKP) K-POINT
C INTERSECTS AT LEAST ONE OF THE EXISTING TETRAHEDRONS:
C WE TRY THE NEXT K-POINT. (GOTO 2000)
C
 3100 continue
      ik=nkadd(nkp)
      volnew=vt
      goto 2500
C
C HERE, WE GET THE NEW TETRAHEDRON WHICH IS MADE OF J-TH FACE
C AND IK-TH K-POINT.
C
 2000 continue
      goto 1100
C
C HERE, WE COULD NOT FIND ANY NEW TETRAHEDRON USING J-TH FACE.
C WE TRY THE NEXT FACE. (GOTO 1100)
C
 2500 continue
      nt=nt+1
      do 2300 i=1,4
      ntetra(i,nt)=ntetra(i,it)
 2300 continue
      ntetra(j,nt)=ik
      voltet(nt)=abs(volnew)
 1100 continue
      if ( it.lt.nt ) goto 1000
      IF (nt>ndiv3 )  CALL juDFT_error(" nt>ndiv3",calledby ="tetcon")
C
C DETERMINE THE CHARACTERISTICS OF THIS DIVISION INTO
C TETRAHEDRA .
C
      vav=0.0
      vsq=0.0
      sav=0.0
      ssq=0.0
      vmin=10.0**9
      vmax=-vmin
      smin=vmin
      smax=vmax
      do it=1,nt
        vt=voltet(it)
        vav=vav+vt
        vsq=vsq+vt**2
        if ( vt.gt.vmax ) vmax=vt
        if ( vt.lt.vmin ) vmin=vt
        do i=1,3
          do j=i+1,4
            dist=0.0
            do icom=1,3
              dist=dist+(kvc(icom,ntetra(i,it))-
     $             kvc(icom,ntetra(j,it)))**2
            end do
            dist=sqrt(dist)
            sav=sav+dist
            ssq=ssq+dist*dist
            if ( dist.gt.smax ) smax=dist
            if ( dist.lt.smin ) smin=dist
          end do
        end do
      end do
      vav=vav/nt
      sav=sav/(6*nt)
      vsq=vsq/nt
      ssq=ssq/(6*nt)
      vsq=sqrt(vsq-vav*vav)
      ssq=sqrt(ssq-sav*sav)
      write(iofile,5000) nt,vav,vsq,vmin,vmax,sav,ssq,smin,smax
 5000 format(/,'   division into tetrahedra  ',/,
     $  '  there are      ',i5,'  tetrahedra ',/,
     $  '  volume         ',f15.10,'  +/-  ',f10.5,3x,2f10.5,/,
     $  '  side           ',f15.10,'  +/-  ',f10.5,3x,2f10.5,/)
c     write(ibfile,5000) nt,vav,vsq,vmin,vmax,sav,ssq,smin,smax
      write(iofile,5100) ((ntetra(j,i),j=1,4),i=1,nt)
 5100 format(4(4x,4i4))
c     write(ibfile,5100) ((ntetra(j,i),j=1,4),i=1,nt)
C CHECK IF WE HAVE THE CORRECT TOTAL VOLUME
      vt=omega*vav*nt*nsym/(2*pi)**3-1.0
      write(iofile,5200) vt
c     write(ibfile,5200) vt
      do 5300 i =1,nt
c     write(ibfile,'('' tetrahedra # '',i5,'' is '',d12.4)') i,voltet(i)
 5300 continue
 5200 format(/,'  voltetsum/volBZ - 1  ',d12.4)
C
C     The following statement used to have a stop in it.
C     If the word TETCON5 appears you have failed the < 1.0D-5 test.
C
      if ( abs(vt).gt.eps1 ) THEN
         write(iofile,'(''  tetcon5  '')')
      endif
      RETURN
      END SUBROUTINE tetcon
      END MODULE m_tetcon

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine trnvec                       *
c     *                                                              *
c     *                       written by : bh                        *
c     *                                                              *
c     *                   last modified : 07/01/90                   *
c     *                                                              *
c     *     this subroutine transforms element vectors to and from   *
c     *     constraint compatable global coordinates and uniform     *
c     *     global coordinates. the element vectors are for one      *
c     *     element in a block of similar, non-conflicting elements. *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine trnvec(vec,trnmte,trne,ndof,nnode,bele,switch)
      implicit integer (a-z)
      include 'param_def'
      double precision
     & vec(mxvl,*), trnmte(mxvl,mxedof,*), sums(mxndof)
      logical trne(mxvl,*)
c
c
c                       note: ndof is an element dependent variable,
c                       but since ndof= 3 for every element in the
c                       program, this is implicitly assumed in this
c                       code. if ever a conflicting element is added
c                       to the prgram, the code will need to be mod-
c                       ified.
c
c
c                       if the transformation matrix for a node is the
c                       identity matrix, skip the transformation. the
c                       the values of switch and type indicate the
c                       following:
c
c                          switch= 1: transformation from uniform global
c                                     coords. to constraint compatable
c                                     global coords.
c
c                          switch= 2: transformation from constraint
c                                     compatable global coords. to
c                                     uniform global coords.
c
c                       the transformation matrices are defined as
c                       taking displacements in uniform global coords.
c                       to those in constraint compatable global coords.
c                       if one transforms from constraint compatable
c                       global coords. to uniform global coords., i.e.
c                       the reverse transformation, one uses the
c                       transpose of the matrices used in the original
c                       transformation. thus the following matrices are
c                       used for the two values of switch:
c
c                          switch= 1 : trans. matrices as stored
c
c                          switch= 2 : transpose
c
      do nod = 1,nnode
         if( trne(bele,nod) ) then
            ts = (nod-1)*ndof
            if( switch.eq.1 ) then
               sums(1)= trnmte(bele,ts+1,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+1,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+1,3)*vec(bele,2*nnode+nod)
               sums(2)= trnmte(bele,ts+2,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+2,3)*vec(bele,2*nnode+nod)
               sums(3)= trnmte(bele,ts+3,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+3,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,3)*vec(bele,2*nnode+nod)
            else
               sums(1)= trnmte(bele,ts+1,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,1)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,1)*vec(bele,2*nnode+nod)
               sums(2)= trnmte(bele,ts+1,2)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,2)*vec(bele,2*nnode+nod)
               sums(3)= trnmte(bele,ts+1,3)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,3)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,3)*vec(bele,2*nnode+nod)
            end if
            vec(bele,nod)         = sums(1)
            vec(bele,nnode+nod)   = sums(2)
            vec(bele,2*nnode+nod) = sums(3)
         end if
      end do
c
c
      return
      end
c     ****************************************************************
c     *                                                              *
c     *                      subroutine trnvec                       *
c     *                                                              *
c     *                       written by : bh                        *
c     *                                                              *
c     *                   last modified : 07/01/90                   *
c     *                                                              *
c     *     this subroutine transforms element vectors to and from   *
c     *     constraint compatable global coordinates and uniform     *
c     *     global coordinates. the element vectors are for one      *
c     *     element in a block of similar, non-conflicting elements. *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine trnvecs(vec,trnmte,trne,ndof,nnode,bele,switch,
     &                   nrow_vec)
      implicit integer (a-z)
      include 'param_def'
      double precision
     & vec(nrow_vec,*), trnmte(mxvl,mxedof,*), sums(mxndof)
      logical trne(mxvl,*)
c
c
c                       note: ndof is an element dependent variable,
c                       but since ndof= 3 for every element in the
c                       program, this is implicitly assumed in this
c                       code. if ever a conflicting element is added
c                       to the prgram, the code will need to be mod-
c                       ified.
c
c
c                       if the transformation matrix for a node is the
c                       identity matrix, skip the transformation. the
c                       the values of switch and type indicate the
c                       following:
c
c                          switch= 1: transformation from uniform global
c                                     coords. to constraint compatable
c                                     global coords.
c
c                          switch= 2: transformation from constraint
c                                     compatable global coords. to
c                                     uniform global coords.
c
c                       the transformation matrices are defined as
c                       taking displacements in uniform global coords.
c                       to those in constraint compatable global coords.
c                       if one transforms from constraint compatable
c                       global coords. to uniform global coords., i.e.
c                       the reverse transformation, one uses the
c                       transpose of the matrices used in the original
c                       transformation. thus the following matrices are
c                       used for the two values of switch:
c
c                          switch= 1 : trans. matrices as stored
c
c                          switch= 2 : transpose
c
      do nod = 1,nnode
         if( trne(bele,nod) ) then
            ts = (nod-1)*ndof
            if( switch.eq.1 ) then
               sums(1)= trnmte(bele,ts+1,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+1,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+1,3)*vec(bele,2*nnode+nod)
               sums(2)= trnmte(bele,ts+2,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+2,3)*vec(bele,2*nnode+nod)
               sums(3)= trnmte(bele,ts+3,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+3,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,3)*vec(bele,2*nnode+nod)
            else
               sums(1)= trnmte(bele,ts+1,1)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,1)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,1)*vec(bele,2*nnode+nod)
               sums(2)= trnmte(bele,ts+1,2)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,2)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,2)*vec(bele,2*nnode+nod)
               sums(3)= trnmte(bele,ts+1,3)*vec(bele,nod)+
     &                  trnmte(bele,ts+2,3)*vec(bele,nnode+nod)+
     &                  trnmte(bele,ts+3,3)*vec(bele,2*nnode+nod)
            end if
            vec(bele,nod)         = sums(1)
            vec(bele,nnode+nod)   = sums(2)
            vec(bele,2*nnode+nod) = sums(3)
         end if
      end do
c
c
      return
      end

ext: f


text:
c     program DRSBLAS4
c>> 1994-10-19 DRSBLAS4 Krogh  Changes to use M77CON
c>> 1992-05-12 DRSBLAS4 CLL Improved consistency of various versions.
c>> 1992-04-15 CLL
c>> 1992-03-16 CLL
c>> 1991-12-02 DRSBLAS4 CLL
c>> 1991-07-31 DRSBLAS4 CLL
c
c     Demonstrates the use of BLAS subroutines SROTMG, SROTM, SAXPY,
c     and SCOPY to implement an algorithm for solving a linear
c     least squares problem using sequential accumulation of the
c     data and Modified (Fast) Givens orthogonal transformations.
c     YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X)
c
c     ------------------------------------------------------------------
c--S replaces "?": DR?BLAS4, ?AXPY, ?ROTMG, ?ROTM, ?COPY, ?VECP
c     ------------------------------------------------------------------
      integer MC, MC1, MC2, MXY
      parameter ( MC=3, MC1=MC+1, MC2=MC+2, MXY=11 )
      real             XTAB(MXY), YTAB(MXY), W(MC2), RG(MC1,MC2)
      real             C(MC), ZERO(1), X, Y, ESTSD, DIV, PARAM(5)
      real             ONE(1)
      integer NXY, NC, NC1, NC2, IXY, J
c
      data XTAB / 0.0e0, .1e0, .2e0, .3e0, .4e0, .5e0, .6e0, .7e0,
     *           .8e0, .9e0, 1.0e0 /
      data YTAB / 1.00e0, .91e0, .86e0, .82e0, .81e0, .82e0, .85e0,
     * .89e0, .95e0, 1.02e0, 1.10e0 /
      data NXY, NC / 11, 3 /
      data ONE(1), ZERO(1) / 1.0e0, 0.0e0 /
c     ------------------------------------------------------------------
      print'(a/)',' DRSBLAS4..  Demo for SROTMG and SROTM.'
      NC1 = NC + 1
      NC2 = NC1 + 1
      call SCOPY(MC1*MC2, ZERO, 0, RG, 1)
      call SCOPY(NC1, ONE, 0, RG(1,NC2), 1)
      do 20 IXY = 1, NXY
        X = XTAB(IXY)
        Y = YTAB(IXY)
c
c          Set W() to be the next row of the least-squares problem.
c
        W(1) = 1.0e0
        W(2) = X
        W(3) = EXP(-X)
        W(4) = Y
        W(5) = 1.0E0
c
c        Accumulate W() into the triangular matrix [R:G]
c        using Givens rotations.
c
        do 10 J = 1, NC1
          call SROTMG(RG(J,NC2), W(NC2), RG(J,J), W(J), PARAM)
          if(J .lt. NC1) call SROTM(NC1-J,RG(J,J+1),MC1,W(J+1),1,PARAM)
   10   continue
c
   20 continue
c
c               Solve the triangular system.
c
        call SCOPY(NC,RG(1,NC1),1,C,1)
        do 30 J = NC, 1, -1
          DIV = RG(J,J)
          if (DIV .EQ. 0.0e0) then
            print*,'ERROR:ZERO DIVISOR AT J =',J
            stop
          end if
          C(J) = C(J) / DIV
          call SAXPY(J-1,-C(J),RG(1,J),1,C,1)
   30 continue
c
      call SVECP(C,NC,'0 Solution: C()=')
      print*,' '
      ESTSD = abs(RG(NC1,NC1)) * sqrt(RG(NC1,NC2)) / sqrt(REAL(NXY-NC))
      print*,'Estimated STD. DEV. of data errors =',ESTSD
      stop
      end


text:
!---------------------------------------------------------------------!
!                                PHAML                                !
!                                                                     !
! The Parallel Hierarchical Adaptive MultiLevel code for solving      !
! linear elliptic partial differential equations of the form          !
! (PUx)x + (QUy)y + RU = F on 2D polygonal domains with mixed         !
! boundary conditions, and eigenvalue problems where F is lambda*U.   !
!                                                                     !
! PHAML is public domain software.  It was produced as part of work   !
! done by the U.S. Government, and is not subject to copyright in     !
! the United States.                                                  !
!                                                                     !
!     William F. Mitchell                                             !
!     Applied and Computational Mathematics Division                  !
!     National Institute of Standards and Technology                  !
!     william.mitchell@nist.gov                                       !
!     http://math.nist.gov/phaml                                      !
!                                                                     !
!---------------------------------------------------------------------!

module phaml_type_mod

!----------------------------------------------------
! This module contains the definition of the primary type, phaml_solution_type.
! Although public, the internals should be considered to be private except
! in module phaml.
!----------------------------------------------------

use global
use message_passing
use gridtype_mod
use zoltan_interf

implicit none
private
public phaml_solution_type

!----------------------------------------------------
! Types defined are :

type phaml_solution_type
   type (grid_type) :: grid
   type (proc_info) :: procs
   integer :: outunit, errunit, pde_id, system_size, eq_type
   character(len=HOSTLEN) :: graphics_host
   type(Zoltan_Struct), pointer :: lb
   character(len=FN_LEN) :: zoltan_param_file
   logical :: i_draw_grid, master_draws_grid,still_sequential
end type phaml_solution_type

end module phaml_type_mod

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
module tensor_rank3_test_module
    use, intrinsic :: iso_fortran_env, only : int64
    use :: util_api, only : assert
    use :: tensor_rank3_module, only : tensor_rank3

    use :: memory_api, only : memory_factory
    use :: tensor_datatypes, only : dt_r64
    use :: tensor_factory_module, only : &
            get_tensor_from_components, &
            create_tensor_from_components
    use :: tensor_module, only : tensor

    implicit none
    private

    public :: tensor_rank3_test

    type :: tensor_rank3_test
    contains
        procedure :: run => run
        procedure :: cleanup => cleanup
        procedure :: clear => clear
    end type tensor_rank3_test

    interface tensor_rank3_test
        module procedure constructor
    end interface tensor_rank3_test
contains
    function constructor() result(this)
        type(tensor_rank3_test) :: this

        call this%clear()
    end function constructor

    subroutine run(this, assertion)
        class(tensor_rank3_test), intent(in) :: this
        type(assert), intent(inout) :: assertion

        type(tensor_rank3) :: atensor_rank3
        class(tensor), allocatable :: atensor

        call assertion%equal("tensor_rank3::Test complete", .true.)

        atensor_rank3 = tensor_rank3(memory_factory%get_data_storage(240), dt_r64, 5, 6,7)

        call assertion%equal("tensor_rank3::Allocated storage", allocated(atensor_rank3%storage))
        call assertion%equal("tensor_rank3::Correct number of bytes", &
                atensor_rank3%storage%number_of_bytes==240)
        call assertion%equal("tensor_rank3::Correct number of elements", &
                atensor_rank3%get_number_of_elements()==210)
        call assertion%equal("tensor_rank3::Correct dimensions", &
                all(atensor_rank3%get_dims()==[5_int64, 6_int64, 7_int64]))
        call assertion%equal("tensor_rank3::Correct rank", &
                atensor_rank3%get_rank()==3)
        call assertion%equal("tensor_rank3::Correct datatype", &
                atensor_rank3%get_datatype()==dt_r64)
        call atensor_rank3%cleanup()

        call create_tensor_from_components(atensor_rank3, memory_factory%get_data_storage(240), dt_r64, &
                [5_int64, 6_int64, 7_int64])

        call assertion%equal("tensor_rank3_from_create_tensor_from_components::Allocated storage", &
                allocated(atensor_rank3%storage))
        call assertion%equal("tensor_rank3_from_create_tensor_from_components::Correct number of bytes", &
                atensor_rank3%storage%number_of_bytes==240)
        call assertion%equal("tensor_rank3_from_create_tensor_from_components::Correct number of elements", &
                atensor_rank3%get_number_of_elements()==210)
        call assertion%equal("tensor_rank3_from_create_tensor_from_components::Correct dimensions", &
                all(atensor_rank3%get_dims()==[5_int64, 6_int64, 7_int64]))
        call assertion%equal("tensor_rank3_from_create_tensor_from_components::Correct rank", &
                atensor_rank3%get_rank()==3)
        call assertion%equal("tensor_rank3_from_create_tensor_from_components::Correct datatype", &
                atensor_rank3%get_datatype()==dt_r64)
        call atensor_rank3%cleanup()

        atensor = get_tensor_from_components(memory_factory%get_data_storage(240), dt_r64, [5_int64, 6_int64, 7_int64])
        call assertion%equal("tensor_rank3_from_get_tensor_from_components::Allocated storage", allocated(atensor%storage))
        call assertion%equal("tensor_rank3_from_get_tensor_from_components::Correct number of bytes", &
                atensor%storage%number_of_bytes==240)
        call assertion%equal("tensor_rank3_from_get_tensor_from_components::Correct number of elements", &
                atensor%get_number_of_elements()==210)
        call assertion%equal("tensor_rank3_from_get_tensor_from_components::Correct dimensions", &
                all(atensor%get_dims()==[5_int64, 6_int64, 7_int64]))
        call assertion%equal("tensor_rank3_from_get_tensor_from_components::Correct rank", &
                atensor%get_rank()==3)
        call assertion%equal("tensor_rank3_from_get_tensor_from_components::Correct datatype", &
                atensor%get_datatype()==dt_r64)
        call atensor%cleanup()
        deallocate(atensor)
    end subroutine run

    subroutine cleanup(this)
        class(tensor_rank3_test), intent(inout) :: this

        call this%clear()
    end subroutine cleanup

    subroutine clear(this)
        class(tensor_rank3_test), intent(inout) :: this
    end subroutine clear
end module tensor_rank3_test_module

max_stars_repo_name: cheshyre/ntcl-tensor

max_issues_repo_name: cheshyre/ntcl-tensor

max_forks_repo_name: cheshyre/ntcl-tensor


text:
      SUBROUTINE TestWout(rzl_array, br, bz, bsupu, bsupv)
      USE stel_kinds, ONLY: rprec, dp
      USE stel_constants, ONLY: twopi
      USE read_wout_mod, ONLY: read_wout_file, read_wout_deallocate,
     1                         tosuvspace, currumnc, currvmnc,
     2                         jcuru, jcurv, xm_nyq, xn_nyq,
     3                         ns_w => ns, mnmax_w => mnmax_nyq,
     4                         xm_w => xm_nyq, xn_w => xn_nyq
      USE vmec_params, ONLY: mscale, nscale, ntmax
      USE vmec_dim, ONLY: ns, ntheta3, ntheta1, nrzt
      USE vmec_input, ONLY: mpol, nzeta, ntor, lasym, nfp
      USE vmec_main, ONLY: hs, lthreed
      USE realspace, ONLY: r1, z1, shalf, sqrts
      USE vforces
      USE vmec_utils, ONLY: getbcyl
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec), INTENT(in)  ::
     1             rzl_array(ns,0:ntor,0:mpol,2*ntmax)
      REAL(rprec), DIMENSION(ns,nzeta,ntheta3), INTENT(in)
     1                         :: bsupu, bsupv
      REAL(rprec), DIMENSION(nrzt), INTENT(in) :: br, bz
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: j, js, k, l, istat, mn, mn0
      REAL(rprec) :: si, ui, vi, r12, z12,
     1               sflx, uflx, rfull, zfull
      REAL(rprec) :: br12, bz12, bphi12,
     1               br1f, bz1f, bphi1f, br2f, bz2f, bphi2f
      LOGICAL :: ls_mesh = .false.               !=T, full mesh, =F, half mesh
C-----------------------------------------------
!      CALL read_wout_file("tichman", istat)
      CALL read_wout_file("fix", istat)
!      CALL read_wout_file("atf14", istat)

!
!     CHECK ON currumn, currvmn vs jcuru, jcurv
!
      DO mn = 1, mnmax_w
         IF (NINT(xm_w(mn)).eq.0 .and. NINT(xn_w(mn)).eq.0) THEN
            mn0 = mn
            EXIT
         END IF
      END DO

      WRITE (30, *)' mnmax_nyq = ', mnmax_w,' ns = ', ns
      WRITE (30, *)
      WRITE (30, *)'  JS     JCURU       CURRU       JCURV       CURRV'
      DO js = 2, ns_w
         WRITE (30, 123) js, jcuru(js), currumnc(mn0,js),
     1                       jcurv(js), currvmnc(mn0,js)
      END DO
 123  FORMAT (i5, 1p,4e12.3)

      RETURN
      IF (istat .ne. 0) RETURN

!      gsqrt => azmn_o    !gsqrt
!      r12 => armn_o; bphi => czmn_o
!      bsupu => crmn_e
!      bsupv => czmn_e

      js = 3*ns/4
!      js = 2                   !!Check this: not working as well as it should at js=2, js>=3 is EXCELLENT!!!
      si = hs*(js-1.5_dp)

      IF (ls_mesh) THEN
        WRITE(33,'(3a)')
     1  '   L  sflx-del  uflx-del      br1/2     bphi1/2       bz1/2',
     2  '     brf-wout  bphif-wout    bzf-wout    brf-vmec  bphif-vmec',
     3  '    bzf-vmec'
      END IF

      DO k = 1, nzeta
         vi = (twopi*(k-1))/nzeta
         DO j = 1, ntheta3
            ui = (twopi*(j-1))/ntheta1
            l = js + ns*(k-1 + nzeta*(j-1))

            IF (ls_mesh) THEN

               rfull = r1(l,0) + sqrts(l)*r1(l,1)
               zfull = z1(l,0) + sqrts(l)*z1(l,1)
!
!          TEST WOUT-BASED VERSION
!
               CALL GetBcyl(rfull, vi/nfp, zfull, br1f, bphi1f, bz1f,
     1                      sflx, uflx, istat)
               IF (istat .ne. 0) PRINT *,' WOUT-BASED ISTAT = ', istat
!
!          TEST VMEC-BASED VERSION
!
               CALL GetBcyl(rfull, vi/nfp, zfull, br2f, bphi2f, bz2f,
     1           sflx, uflx, bsupu, bsupv, rzl_array, ns, ntor, mpol,
     2           ntmax, nzeta, ntheta3, nfp, mscale, nscale,
     3           lthreed, lasym, istat)
               IF (istat .ne. 0) PRINT *,' VMEC-BASED ISTAT = ', istat

               WRITE (33, 1225) l, sflx-si, MOD(uflx-ui,twopi),
     1            br(l), czmn_o(l), bz(l),
     2            br1f, bphi1f, bz1f, br2f, bphi2f, bz1f

 1225   FORMAT(i4, 1p,2e10.2, 9e12.3)

            ELSE

            z12 = (z1(l,0) + z1(l-1,0) + shalf(l)*
     1            (z1(l,1) + z1(l-1,1)))/2
            r12 = armn_o(l)

!
!           TEST WOUT-BASED VERSION
!
            CALL GetBcyl(r12, vi/nfp, z12, br12, bphi12, bz12,
     1                   sflx, uflx, istat)

            IF (istat .eq. 0) THEN
               WRITE (33, 1224) l, sflx-si, MOD(uflx-ui,twopi),
     1            (br(l) - br12)/bphi12, (czmn_o(l)-bphi12)/bphi12,
     2            (bz(l) - bz12)/bphi12
            ELSE
               WRITE (33, *)' istat = ', istat
            END IF

!
!          TEST VMEC-BASED VERSION (OVERLOADED, SAME FUNCTION CALL, DIFF ARG.)
!
            CALL GetBcyl(r12, vi/nfp, z12, br12, bphi12, bz12, sflx,
     1           uflx, bsupu, bsupv, rzl_array, ns, ntor, mpol,
     2           ntmax, nzeta, ntheta3, nfp, mscale, nscale,
     3           lthreed, lasym, istat)
            IF (istat .eq. 0) THEN
               WRITE (34, 1224) l, sflx-si, MOD(uflx-ui,twopi),
     1            (br(l) - br12)/bphi12, (czmn_o(l)-bphi12)/bphi12,
     2            (bz(l) - bz12)/bphi12
            ELSE
               WRITE (34, *)' istat = ', istat
            END IF

            END IF

!            CALL tosuvspace(si,ui,vi,gsqrt=gsqrt1)
!            WRITE(33, 1223) l, gsqrt(l), gsqrt1,
!     1                 ABS(gsqrt(l)-gsqrt1)/ABS(gsqrt1)
         END DO
      END DO
 1223 FORMAT(' l = ',i4,' GSQRT = ',1p,e12.4,' FROM WOUT = ',e12.4,
     1      ' REL.ERROR = ',e10.2)
 1224 FORMAT(' l = ',i4, ' del-SFLX = ', 1p,e10.2,
     1                   ' del-UFLX = ', e10.2,
     2                   ' del-BR = ', e12.4,' del-BPHI = ',
     3            e12.4, ' del-BZ = ', e12.4)

      CALL read_wout_deallocate

      END SUBROUTINE TestWout

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT


text:
      REAL*8 FUNCTION RPMAX2(R1,R2,M1,M2,KS1,KS2,VINF)
*
*
*       Maximum periastron factor for capture.
*       --------------------------------------
*
*       Developed by Rosemary Mardling (March 1996).
*       @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*
      REAL*8  R1,R2,M1,M2,VINF,M21,A(8),B(2),P(2)
      INTEGER KS1,KS2,KS(2)
      DATA  A /-600.2227,3538.748,1703.075,77.31540,
     &         191.8260,-2.203482,1202.116,70.80281/
      DATA B /0.54,1.08/
*
*
*       Define mass and radius ratios and swap stellar type if necessary.
      IF (R1.GE.R2) THEN
          M21 = M2/M1
          R21 = R2/R1
          KS(1) = KS1
	  KS(2) = KS2
      ELSE
	  M21 = M1/M2
	  R21 = R1/R2
	  KS(1) = KS2
	  KS(2) = KS1
      END IF
*
*       Choose polytropic index for each star (adopt n=1.5 for giants).
      DO 1 I = 1,2
          IF(KS(I).EQ.0.OR.KS(I).EQ.3.OR.KS(I).EQ.5.OR.
     &       KS(I).EQ.6.OR.KS(I).EQ.9)THEN
              IP = 1
              PP = B(1)
          ELSE
              IP = 5
              PP = B(2)
          ENDIF
          P(I) = (((A(IP+3)*M21 + A(IP+2))*M21 + A(IP+1))*M21
     &                                         + A(IP))/VINF**PP
    1 CONTINUE
*
*       Calculate maximum periastron separation for capture (note arg < 0).
      Z = P(1) + (R21**5/M21**2)*P(2)
      IF (Z.LE.0.0) THEN
*       Adopt factor 2 for collision test in case of small mass ratio.
          RPMAX2 = 2.0
      ELSE
          RPMAX2 = Z**(1.0/6.0)
*         RPMAX2 = (P(1) + (R21**5/M21**2)*P(2))**(1.0/6.0)
      END IF
*
      RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!===============================================================================
! One of Andy Nowacki's Fortran utility modules for dealing with seismic
! anisotropy and other problems.
!
! Andy Nowacki <andy.nowacki@bristol.ac.uk>
!
! See the file LICENCE for licence details.
!===============================================================================
   module density_1d
!------------------------------------------------------------------------------
!
! Returns density for some global Earth models, given a depth.
!
!
!  ** size constants
      integer, parameter, private :: i4 = selected_int_kind(9) ; ! long int
      integer, parameter, private :: r4 = selected_real_kind(6,37) ; ! SP
      integer, parameter, private :: r8 = selected_real_kind(15,307) ; ! DP

!  ** precision selector
      integer, parameter, private :: rs = r8

!  ** maths constants and other useful things
      real(rs), parameter, private :: pi = 3.141592653589793238462643 ;
      real(rs), parameter, private :: to_rad = 1.74532925199433e-002 ;
      real(rs), parameter, private :: to_deg = 57.2957795130823e0 ;
      real(rs), parameter, private :: to_km = 111.194926644559 ;

      real(rs), parameter, private :: big_number = 10.e36 ;

   contains

!===============================================================================
   function ak135den(depth)
!===============================================================================
!  Uses AK135 model

   implicit none

   real(rs),intent(in)  :: depth
   real(rs)  :: dep(136),den(136),ak135den,den1,den2,dep1,dep2
   integer   :: i
   logical   :: twodepths
   data dep / &
		 0.0,   20.0,   20.0,   35.0,   35.0, &
		77.5,  120.0,  165.0,  210.0,  210.0, &
	   260.0,  310.0,  360.0,  410.0,  410.0, &
	   460.0,  510.0,  560.0,  610.0,  660.0, &
	   660.0,  710.0,  760.0,  809.5,  859.0, &
	   908.5,  958.0, 1007.5, 1057.0, 1106.5, &
	  1156.0, 1205.5, 1255.0, 1304.5, 1354.0, &
	  1403.5, 1453.0, 1502.5, 1552.0, 1601.5, &
	  1651.0, 1700.5, 1750.0, 1799.5, 1849.0, &
	  1898.5, 1948.0, 1997.5, 2047.0, 2096.5, &
	  2146.0, 2195.5, 2245.0, 2294.5, 2344.0, &
	  2393.5, 2443.0, 2492.5, 2542.0, 2591.5, &
	  2640.0, 2690.0, 2740.0, 2740.0, 2789.7, &
	  2839.3, 2891.5, 2891.5, 2939.3, 2989.7, &
	  3040.0, 3090.3, 3140.7, 3191.0, 3241.3, &
	  3291.7, 3342.0, 3392.3, 3442.6, 3493.0, &
	  3543.3, 3593.6, 3644.0, 3694.3, 3744.6, &
	  3795.0, 3845.3, 3895.6, 3945.9, 3996.3, &
	  4046.6, 4096.9, 4147.3, 4197.6, 4247.9, &
	  4298.3, 4348.6, 4398.9, 4449.3, 4499.6, &
	  4549.9, 4600.3, 4650.6, 4700.9, 4801.6, &
	  4851.9, 4902.2, 4952.6, 5002.9, 5053.2, &
	  5103.6, 5153.5, 5153.5, 5204.6, 5255.3, &
	  5306.0, 5356.8, 5407.5, 5458.2, 5508.9, &
	  5559.6, 5610.3, 5661.0, 5711.7, 5813.2, &
	  5863.9, 5914.6, 5965.3, 6016.0, 6066.7, &
	  6117.4, 6168.1, 6218.9, 6269.6, 6320.3, &
	  6371.0 /

   data den / &
	   2.44900,  2.44900,  2.71420,  2.71420,  3.29760, &
	   3.29940,  3.30130,  3.34870,  3.39600,  3.39600, &
	   3.46520,  3.53430,  3.60340,  3.67260,  3.79760, &
	   3.86120,  3.92480,  3.98850,  4.05210,  4.11580, &
	   4.33930,  4.38960,  4.43990,  4.47010,  4.50290, &
	   4.53500,  4.56640,  4.59700,  4.62700,  4.65630, &
	   4.68490,  4.71300,  4.74040,  4.76730,  4.79360, &
	   4.81950,  4.84480,  4.87000,  4.89420,  4.91820, &
	   4.94230,  4.96530,  4.98820,  5.01090,  5.03330, &
	   5.05530,  5.07720,  5.09900,  5.12060,  5.14220, &
	   5.16380,  5.18480,  5.20610,  5.22700,  5.24810, &
	   5.26980,  5.29070,  5.31220,  5.33380,  5.35600, &
	   5.37760,  5.39990,  5.42240,  5.42240,  5.42380, &
	   5.42510,  5.42650,  9.91450,  9.99420, 10.07220, &
	  10.14850, 10.22330, 10.29640, 10.36790, 10.43780, &
	  10.50620, 10.57310, 10.63850, 10.70230, 10.76470, &
	  10.82570, 10.88520, 10.94340, 11.00010, 11.05550, &
	  11.10950, 11.16230, 11.21370, 11.26390, 11.31270, &
	  11.36040, 11.40690, 11.45210, 11.49620, 11.53910, &
	  11.58090, 11.62160, 11.66120, 11.69980, 11.73730, &
	  11.77370, 11.80920, 11.84370, 11.87720, 11.94140, &
	  11.97220, 12.00010, 12.03110, 12.05930, 12.08670, &
	  12.11330, 12.13910, 12.70370, 12.72890, 12.75300, &
	  12.77600, 12.79800, 12.81880, 12.83870, 12.85740, &
	  12.87510, 12.89170, 12.90720, 12.92170, 12.94740, &
	  12.95860, 12.96880, 12.97790, 12.98590, 12.99290, &
	  12.99880, 13.00360, 13.00740, 13.01000, 13.01170, &
	  13.01220 /

!  Ensure depth is in km and not too large
   if (depth > 6371.0) then
      write(*,'(a,f9.1,a)') '  ak135den: depth ',depth,' is too large.'
      stop
   endif

!  Look for the specific point first
   do i=1,136
      if (depth == dep(i)) then
         ak135den = den(i)
         twodepths = .false.
         exit
      else
!  Otherwise, find the points above and below this depth
         if (depth > dep(i) .and. depth < dep(i+1)) then
            dep1 = dep(i)
            den1 = den(i)
            dep2 = dep(i+1)
            den2 = den(i+1)
            twodepths = .true.
            exit
         endif
      endif
   enddo

   if (twodepths) then
!  Extrapolate linearly to find density
      ak135den = den1 + (depth-dep1)*(den2-den1)/(dep2-dep1)
   endif

   return

   end function ak135den
!------------------------------------------------------------------------------

!===============================================================================
   function premden(depth)
!===============================================================================
!  PREM

   implicit none

   real(rs),intent(in) :: depth
   real(rs)  :: dep(88),den(88),premden,den1,den2,dep1,dep2
   integer   :: i
   logical   :: twodepths
   data dep /        0.00,   15.00,   15.00,   24.40,   24.40, &
		   40.00,   60.00,   80.00,  115.00,  150.00,  185.00, &
		  220.00,  220.00,  265.00,  310.00,  355.00,  400.00, &
		  400.00,  450.00,  500.00,  550.00,  600.00,  635.00, &
		  670.00,  670.00,  721.00,  771.00,  871.00,  971.00, &
		 1071.00, 1171.00, 1271.00, 1371.00, 1471.00, 1571.00, &
		 1671.00, 1771.00, 1871.00, 1971.00, 2071.00, 2171.00, &
		 2271.00, 2371.00, 2471.00, 2571.00, 2671.00, 2741.00, &
		 2771.00, 2871.00, 2891.00, 2891.00, 2971.00, 3071.00, &
		 3171.00, 3271.00, 3371.00, 3471.00, 3571.00, 3671.00, &
		 3771.00, 3871.00, 3971.00, 4017.00, 4171.00, 4271.00, &
		 4371.00, 4471.00, 4571.00, 4671.00, 4771.00, 4871.00, &
		 4971.00, 5071.00, 5149.50, 5149.50, 5171.00, 5271.00, &
		 5371.00, 5471.00, 5571.00, 5671.00, 5771.00, 5871.00, &
		 5971.00, 6071.00, 6171.00, 6271.00, 6371.00/

   data den /       2.60000,  2.60000,  2.90000,  2.90000,  3.38076, &
		  3.37906,  3.37688,  3.37471,  3.37091,  3.36710,  3.36330, &
		  3.35950,  3.43578,  3.46264,  3.48951,  3.51639,  3.54325, &
		  3.72378,  3.78678,  3.84980,  3.91282,  3.97584,  3.98399, &
		  3.99214,  4.38071,  4.41241,  4.44317,  4.50372,  4.56307, &
		  4.62129,  4.67844,  4.73460,  4.78983,  4.84422,  4.89783, &
		  4.95073,  5.00299,  5.05469,  5.10590,  5.15669,  5.20713, &
		  5.25729,  5.30724,  5.35706,  5.40681,  5.45657,  5.49145, &
		  5.50642,  5.55641,  5.56645,  9.90349, 10.02940, 10.18134, &
		 10.32726, 10.46727, 10.60152, 10.73012, 10.85321, 10.97091, &
		 11.08335, 11.19067, 11.29298, 11.39042, 11.48311, 11.57119, &
		 11.65478, 11.73401, 11.80900, 11.87990, 11.94682, 12.00989, &
		 12.06924, 12.12500, 12.16634, 12.76360, 12.77493, 12.82501, &
		 12.87073, 12.91211, 12.94912, 12.98178, 13.01009, 13.03404, &
		 13.05364, 13.06888, 13.07977, 13.08630, 13.08848 /

!  Ensure depth is in km and not too large
   if (depth > 6371.0) then
      write(*,'(a,f9.1,a)') '  premden: depth ',depth,' is too large.'
      stop
   endif

!  Look for the specific point first
   do i=1,88
      if (depth == dep(i)) then
         premden = den(i)
         twodepths = .false.
         exit
      else
!  Otherwise, find the points above and below this depth
         if (depth > dep(i) .and. depth < dep(i+1)) then
            dep1 = dep(i)
            den1 = den(i)
            dep2 = dep(i+1)
            den2 = den(i+1)
            twodepths = .true.
            exit
         endif
      endif
   enddo

   if (twodepths) then
!  Extrapolate linearly to find density
      premden = den1 + (depth-dep1)*(den2-den1)/(dep2-dep1)
   endif

   return

   end function premden
!------------------------------------------------------------------------------

   end module density_1d

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork