text:
      subroutine dirmago2op(socsc,e,l,my,vr,bspr,bopr,
     >                    dx,xnot,rs,ns,g,f,gp,iprpts)
c
      implicit real*8(a-h,o-z)
!     include '../param.h'
!     include 'atom_param.h'
c
c     ***********************************************************
c     * integration of relativistic radial dirac equation       *
c     * in the presence of an internal field by adams method    *
c     * integrate outward!                                      *
c     *                                                         *
c     * strange et al., j.phys.c.solid state phys.17,3355(1984) *
c     * scaling of SOC and OP-term ala Hubert Ebert included    *
c     ***********************************************************
c
      integer iprpts

      complex*16 e,p1,q1,p2,q2,pp1,qp1,pp2,qp2
      complex*16 p1n,q1n,k11,k12,k13,k14,m11,m12,m13,m14
      complex*16 p2n,q2n,k21,k22,k23,k24,m21,m22,m23,m24
      complex*16 p1sum,q1sum,p10,q10,p11,q11,pp1np1,qp1np1
      complex*16 p2sum,q2sum,p20,q20,p21,q21,pp2np1,qp2np1
      complex*16 g,gp,f
c
      dimension vr(iprpts),bspr(iprpts),bopr(iprpts,2)
      dimension ba11(iprpts),ba12(iprpts),ba22(iprpts)
      dimension bb11(iprpts),bb22(iprpts)
      dimension p1(iprpts),p2(iprpts),q1(iprpts),q2(iprpts)
      dimension pp1(iprpts),pp2(iprpts),qp1(iprpts),qp2(iprpts)
      dimension fk1(2),fk2(2),gk1(2),gk2(2)
      dimension g(2,2,iprpts),gp(2,2),f(2,2,iprpts)
      data nitmax/100/,test/1.0d10/,imode/0/
c
      if(iabs(my).eq.2*l+1) stop 'dirmag2: my=+/-(l+1/2)'
      dkoef1=475.d0/502.d0
      dkoef2= 27.d0/502.d0
      xl = dfloat(l)
      dx2=0.5d0*dx
      c = 274.072d0
      cin = 1.0d 00/(c*c)
c
c     calculate boundary conditions
c
      if(vr(1).lt.-1.d-3)then
      hoc=-vr(1)/c
      else
c     write(6,'(99d12.4)')(vr(i),i=1,5)
      hoc=-vr(2)/c
      vr(1)=vr(2)
      endif

c
      xl=dfloat(l)
      xkap=xl
      xk1=socsc*(1.d0+xkap)-1.d0
      facl1=xl*(xl+1.d0)-xk1*(xk1+1.d0)
      u1=(xk1+dsqrt(xk1*xk1-hoc*hoc+facl1))/hoc
      xkap=-xl-1.d0
      xk2=socsc*(1.d0+xkap)-1.d0
      facl2=xl*(xl+1.d0)-xk2*(xk2+1.d0)
      u2=(xk2+dsqrt(xk2*xk2-hoc*hoc+facl2))/hoc
c
c     get combined effective fields
c
      call brmat(l,my,ns,bspr,bopr,ba11,ba12,ba22,bb11,bb22,iprpts)
c
      do ii=1,2
      if(ii.eq.2) go to 40
c
      p1(1) = 1.0d-20
      q1(1) = c*u1*1.0d-20
      p2(1) = (0.0d 00,0.0d 00)
      q2(1) = (0.0d 00,0.0d 00)
      go to 50
c
   40 continue
      p1(1) = (0.0d 00,0.0d 00)
      q1(1) = (0.0d 00,0.0d 00)
      p2(1) = 1.0d-20
      q2(1) = c*u2*1.0d-20
   50 continue
c
      r=dexp(xnot)
      call dmag2op(pp1(1),qp1(1),pp2(1),qp2(1),
     >           p1(1),q1(1),p2(1),q2(1),
     >           xk1,xk2,facl1,facl2,cin,e,r,
     >           vr(1),ba11(1),ba12(1),ba22(1),bb11(1),bb22(1))
c
c
c---> start runge-kutta procedure (points 2, ... , 6)
c
      do n=1,5
c
        x=xnot+(n-1)*dx
        rn=dexp(x)
        vrn=vr(n)
        ba11n=ba11(n)
        ba12n=ba12(n)
        ba22n=ba22(n)
        bb11n=bb11(n)
        bb22n=bb22(n)
c
        rnp1=dexp(x+dx)
        vrnp1=vr(n+1)
        ba11np1=ba11(n+1)
        ba12np1=ba12(n+1)
        ba22np1=ba22(n+1)
        bb11np1=bb11(n+1)
        bb22np1=bb22(n+1)
c
        rmid=dexp(x+dx2)
        vrmid=0.5d0*(vrn+vrnp1)
        ba11mid=0.5d0*(ba11n+ba11np1)
        ba12mid=0.5d0*(ba12n+ba12np1)
        ba22mid=0.5d0*(ba22n+ba22np1)
        bb11mid=0.5d0*(bb11n+bb11np1)
        bb22mid=0.5d0*(bb22n+bb22np1)
c
        p1n=p1(n)
        q1n=q1(n)
        p2n=p2(n)
        q2n=q2(n)
        call dmag2op(k11,m11,k21,m21,p1n,q1n,p2n,q2n,
     >             xk1,xk2,facl1,facl2,cin,e,rn,
     >             vrn,ba11n,ba12n,ba22n,bb11n,bb22n)
c
        call dmag2op(k12,m12,k22,m22,
     >             p1n+dx2*k11,q1n+dx2*m11,p2n+dx2*k21,q2n+dx2*m21,
     >             xk1,xk2,facl1,facl2,cin,e,rmid,
     >             vrmid,ba11mid,ba12mid,ba22mid,bb11mid,bb22mid)
c
        call dmag2op(k13,m13,k23,m23,
     >             p1n+dx2*k12,q1n+dx2*m12,p2n+dx2*k22,q2n+dx2*m22,
     >             xk1,xk2,facl1,facl2,cin,e,rmid,
     >             vrmid,ba11mid,ba12mid,ba22mid,bb11mid,bb22mid)
c
        call dmag2op(k14,m14,k24,m24,
     >             p1n+dx*k13,q1n+dx*m13,p2n+dx*k23,q2n+dx*m23,
     >             xk1,xk2,facl1,facl2,cin,e,rnp1,
     >             vrnp1,ba11np1,ba12np1,ba22np1,bb11np1,bb22np1)
c
        p1(n+1)=p1n+dx*(k11+2.d0*k12+2.d0*k13+k14)/6.d0
        q1(n+1)=q1n+dx*(m11+2.d0*m12+2.d0*m13+m14)/6.d0
        p2(n+1)=p2n+dx*(k21+2.d0*k22+2.d0*k23+k24)/6.d0
        q2(n+1)=q2n+dx*(m21+2.d0*m22+2.d0*m23+m24)/6.d0
        call dmag2op(pp1(n+1),qp1(n+1),pp2(n+1),qp2(n+1),
     >             p1(n+1),q1(n+1),p2(n+1),q2(n+1),
     >             xk1,xk2,facl1,facl2,cin,e,rnp1,
     >             vrnp1,ba11np1,ba12np1,ba22np1,bb11np1,bb22np1)
c
      end do
c
c---> begin adams procedure (points 7, 8, ... , ns)
c
      do n=6,ns-1
c
        x=xnot+(n-1)*dx
        r=dexp(x+dx)
        vrc=vr(n+1)
        ba11c=ba11(n+1)
        ba12c=ba12(n+1)
        ba22c=ba22(n+1)
        bb11c=bb11(n+1)
        bb22c=bb22(n+1)
c
        p1sum=646.d0*pp1(n)-264.d0*pp1(n-1)+106.d0*pp1(n-2)-
     >        19.d0*pp1(n-3)
        q1sum=646.d0*qp1(n)-264.d0*qp1(n-1)+106.d0*qp1(n-2)-
     >        19.d0*qp1(n-3)
        p2sum=646.d0*pp2(n)-264.d0*pp2(n-1)+106.d0*pp2(n-2)-
     >        19.d0*pp2(n-3)
        q2sum=646.d0*qp2(n)-264.d0*qp2(n-1)+106.d0*qp2(n-2)-
     >        19.d0*qp2(n-3)
c
c  predict for point n+1
c
        p10=p1(n)+dx*(251.d0*pp1(n-4)-1274.d0*pp1(n-3)+
     >      2616.d0*pp1(n-2)-2774.d0*pp1(n-1)+1901.d0*pp1(n))/720.d0
        q10=q1(n)+dx*(251.d0*qp1(n-4)-1274.d0*qp1(n-3)+
     >      2616.d0*qp1(n-2)-2774.d0*qp1(n-1)+1901.d0*qp1(n))/720.d0
        p20=p2(n)+dx*(251.d0*pp2(n-4)-1274.d0*pp2(n-3)+
     >      2616.d0*pp2(n-2)-2774.d0*pp2(n-1)+1901.d0*pp2(n))/720.d0
        q20=q2(n)+dx*(251.d0*qp2(n-4)-1274.d0*qp2(n-3)+
     >      2616.d0*qp2(n-2)-2774.d0*qp2(n-1)+1901.d0*qp2(n))/720.d0
c
c  correct
c
        if(imode.eq.1) then
c
        do nit=1,nitmax
c
          call dmag2op(pp1np1,qp1np1,pp2np1,qp2np1,
     >               p10,q10,p20,q20,
     >               xk1,xk2,facl1,facl2,cin,e,r,
     >               vrc,ba11c,ba12c,ba22c,bb11c,bb22c)
          p11=p1(n)+dx*(251.d0*pp1np1+p1sum)/720.d0
          q11=q1(n)+dx*(251.d0*qp1np1+q1sum)/720.d0
          p21=p2(n)+dx*(251.d0*pp2np1+p2sum)/720.d0
          q21=q2(n)+dx*(251.d0*qp2np1+q2sum)/720.d0
c
c  compare predictor with corrector
c
          if(test*cdabs(p11-p10).gt.cdabs(p10)) goto 15
          if(test*cdabs(q11-q10).gt.cdabs(q10)) goto 15
          if(test*cdabs(p21-p20).gt.cdabs(p20)) goto 15
          if(test*cdabs(q21-q20).gt.cdabs(q20)) goto 15
          goto 10
c
   15     p10=p11
          q10=q11
          p20=p21
          q20=q21
        end do
        write(6,*) n+1,r,nit,' not converged'
c
        else
c
        call dmag2op(pp1np1,qp1np1,pp2np1,qp2np1,
     >             p10,q10,p20,q20,
     >             xk1,xk2,facl1,facl2,cin,e,r,
     >             vrc,ba11c,ba12c,ba22c,bb11c,bb22c)
        p11=p1(n)+dx*(251.d0*pp1np1+p1sum)/720.d0
        q11=q1(n)+dx*(251.d0*qp1np1+q1sum)/720.d0
        p21=p2(n)+dx*(251.d0*pp2np1+p2sum)/720.d0
        q21=q2(n)+dx*(251.d0*qp2np1+q2sum)/720.d0
c
        q11=dkoef1*q11+dkoef2*q10
        p11=dkoef1*p11+dkoef2*p10
        q21=dkoef1*q21+dkoef2*q20
        p21=dkoef1*p21+dkoef2*p20
c
        end if
c
   10   q1(n+1)=q11
        p1(n+1)=p11
        q2(n+1)=q21
        p2(n+1)=p21
        call dmag2op(pp1(n+1),qp1(n+1),pp2(n+1),qp2(n+1),
     >             p1(n+1),q1(n+1),p2(n+1),q2(n+1),
     >             xk1,xk2,facl1,facl2,cin,e,r,
     >             vrc,ba11c,ba12c,ba22c,bb11c,bb22c)
c
      end do
c
c     store radial amplitudes times radius
c
      do n=1,ns
        g(1,ii,n)=p1(n)
        g(2,ii,n)=p2(n)
        f(1,ii,n)=q1(n)/c
        f(2,ii,n)=q2(n)/c
      end do
      gp(1,ii)=(pp1(ns)-p1(ns))/rs
      gp(2,ii)=(pp2(ns)-p2(ns))/rs
c
      end do
c
      return
      end
      subroutine dirmagi2op(socsc,e,l,my,vr,bspr,bopr,
     >                    dx,xnot,rs,ns,g,f,gp,iprpts)
c
      implicit real*8(a-h,o-z)
!     include '../param.h'
!     include 'atom_param.h'
c
c     ***********************************************************
c     * integration of relativistic radial dirac equation       *
c     * in the presence of an internal field by adams method    *
c     * integrate inward!                                       *
c     *                                                         *
c     * strange et al., j.phys.c.solid state phys.17,3355(1984) *
c     * scaling of SOC and OP-term ala Hubert Ebert included    *
c     ***********************************************************
c
      integer iprpts

      complex*16 e,psq,pe
      complex*16 p1,q1,p2,q2,pp1,qp1,pp2,qp2
      complex*16 p1n,q1n,k11,k12,k13,k14,m11,m12,m13,m14
      complex*16 p2n,q2n,k21,k22,k23,k24,m21,m22,m23,m24
      complex*16 p1sum,q1sum,p10,q10,p11,q11,pp1nm1,qp1nm1
      complex*16 p2sum,q2sum,p20,q20,p21,q21,pp2nm1,qp2nm1
      complex*16 g,gp,f
      complex*16 fb(0:5),fn(0:5),fh(0:5)
c
      dimension vr(iprpts),bspr(iprpts),bopr(iprpts,2)
      dimension xlag(iprpts),vrlag(iprpts)
      dimension bsprlag(iprpts),boprlag(iprpts,2)
      dimension ba11(iprpts),ba12(iprpts)
      dimension ba22(iprpts),bb11(iprpts),bb22(iprpts)
      dimension p1(iprpts),p2(iprpts),q1(iprpts),q2(iprpts)
      dimension pp1(iprpts),pp2(iprpts),qp1(iprpts),qp2(iprpts)
      dimension fk1(2),fk2(2),gk1(2),gk2(2)
      dimension g(2,2,iprpts),gp(2,2),f(2,2,iprpts)
      data nitmax/100/,test/1.0d10/,imode/0/,nout/5/
c
      if(iabs(my).eq.2*l+1) stop 'dirmag2: my=+/-(l+1/2)'
      dkoef1=475.d0/502.d0
      dkoef2= 27.d0/502.d0
c
      lb1=l-1
      lb2=l+1
      sk1=l-lb1
      sk2=l-lb2
c
      xl = dfloat(l)
      xkap=xl
      xk1=socsc*(1.d0+xkap)-1.d0
      facl1=xl*(xl+1.d0)-xk1*(xk1+1.d0)
      xkap=-xl-1.d0
      xk2=socsc*(1.d0+xkap)-1.d0
      facl2=xl*(xl+1.d0)-xk2*(xk2+1.d0)
c
      dx2=0.5d0*dx
      c = 274.072d0
      cin = 1.0d 00/(c*c)
      ilag=3
      xs=dlog(rs)
      do i=1,ns
        xlag(i)=xs+(i-ns)*dx
        vrlag(i)=vr(i)
        bsprlag(i)=bspr(i)
        boprlag(i,1)=bopr(i,1)
        boprlag(i,2)=bopr(i,2)
      end do
      nlag=ns+1
      vrlag(nlag)=0.d0
      bsprlag(nlag)=0.d0
      boprlag(nlag,1)=0.d0
      boprlag(nlag,2)=0.d0
      xlag(nlag)=xs+nout*dx
c
c     get combined effective fields
c
      call brmat(l,my,nlag,bsprlag,boprlag,ba11,ba12,ba22,bb11,bb22,
     >           iprpts)
c
      psq=e+e*e*cin
      pe=cdsqrt(psq)
c
      do ii=1,2
c
c---> set up the starting values outside the muffin tin
c     corresponding to the boundary condition
c
      x=xs+nout*dx
      r=dexp(x)
      call csbf(l+1,pe,r,fb,fn,fh)
      n=ns+nout
c
      if(ii.eq.2) go to 40
      p1(n) = r*fb(l)
      q1(n) = ( (xk1-xl)*fb(l) + sk1*pe*r*fb(lb1) ) * e/psq
      p2(n) = 0.d0
      q2(n) = 0.d0
      call dmag2op(pp1(n),qp1(n),pp2(n),qp2(n),
     >           p1(n),q1(n),p2(n),q2(n),
     >           xk1,xk2,facl1,facl2,cin,e,r,
     >           vrlag(nlag),ba11(nlag),ba12(nlag),ba22(nlag),
     >                       bb11(nlag),bb22(nlag))
      go to 50
   40 continue
      p1(n) = 0.d0
      q1(n) = 0.d0
      p2(n) = r*fb(l)
      q2(n) = ( (xk2+xl+1.d0)*fb(l) + sk2*pe*r*fb(lb2) ) * e/psq
      call dmag2op(pp1(n),qp1(n),pp2(n),qp2(n),
     >           p1(n),q1(n),p2(n),q2(n),
     >           xk1,xk2,facl1,facl2,cin,e,r,
     >           vrlag(nlag),ba11(nlag),ba12(nlag),ba22(nlag),
     >                       bb11(nlag),bb22(nlag))
   50 continue
c
c     write(6,*) n
c     write(6,'(8d15.7)') p1(n),q1(n),p2(n),q2(n)
c     write(6,'(8d15.7)') pp1(n),qp1(n),pp2(n),qp2(n)
c     write(6,*)
c
c---> start runge-kutta procedure (points ns+nout-1, ... , ns+nout-5)
c
      do n=ns+nout,ns+nout-4,-1
c
        x=xs+(n-ns)*dx
        rn=dexp(x)
        vrn=ylag(x,xlag,vrlag,0,ilag,nlag,iex)
        ba11n=ylag(x,xlag,ba11,0,ilag,nlag,iex)
        ba12n=ylag(x,xlag,ba12,0,ilag,nlag,iex)
        ba22n=ylag(x,xlag,ba22,0,ilag,nlag,iex)
        bb11n=ylag(x,xlag,bb11,0,ilag,nlag,iex)
        bb22n=ylag(x,xlag,bb22,0,ilag,nlag,iex)
c
        rmid=dexp(x-dx2)
        vrmid=ylag(x-dx2,xlag,vrlag,0,ilag,nlag,iex)
        ba11mid=ylag(x-dx2,xlag,ba11,0,ilag,nlag,iex)
        ba12mid=ylag(x-dx2,xlag,ba12,0,ilag,nlag,iex)
        ba22mid=ylag(x-dx2,xlag,ba22,0,ilag,nlag,iex)
        bb11mid=ylag(x-dx2,xlag,bb11,0,ilag,nlag,iex)
        bb22mid=ylag(x-dx2,xlag,bb22,0,ilag,nlag,iex)
c
        rnm1=dexp(x-dx)
        vrnm1=ylag(x-dx,xlag,vrlag,0,ilag,nlag,iex)
        ba11nm1=ylag(x-dx,xlag,ba11,0,ilag,nlag,iex)
        ba12nm1=ylag(x-dx,xlag,ba12,0,ilag,nlag,iex)
        ba22nm1=ylag(x-dx,xlag,ba22,0,ilag,nlag,iex)
        bb11nm1=ylag(x-dx,xlag,bb11,0,ilag,nlag,iex)
        bb22nm1=ylag(x-dx,xlag,bb22,0,ilag,nlag,iex)
c
        p1n=p1(n)
        q1n=q1(n)
        p2n=p2(n)
        q2n=q2(n)
        call dmag2op(k11,m11,k21,m21,p1n,q1n,p2n,q2n,
     >             xk1,xk2,facl1,facl2,cin,e,rn,
     >             vrn,ba11n,ba12n,ba22n,bb11n,bb22n)
c
        call dmag2op(k12,m12,k22,m22,
     >             p1n-dx2*k11,q1n-dx2*m11,p2n-dx2*k21,q2n-dx2*m21,
     >             xk1,xk2,facl1,facl2,cin,e,rmid,
     >             vrmid,ba11mid,ba12mid,ba22mid,bb11mid,bb22mid)
c
        call dmag2op(k13,m13,k23,m23,
     >             p1n-dx2*k12,q1n-dx2*m12,p2n-dx2*k22,q2n-dx2*m22,
     >             xk1,xk2,facl1,facl2,cin,e,rmid,
     >             vrmid,ba11mid,ba12mid,ba22mid,bb11mid,bb22mid)
c
        call dmag2op(k14,m14,k24,m24,
     >             p1n-dx*k13,q1n-dx*m13,p2n-dx*k23,q2n-dx*m23,
     >             xk1,xk2,facl1,facl2,cin,e,rnm1,
     >             vrnm1,ba11nm1,ba12nm1,ba22nm1,bb11nm1,bb22nm1)
c
        p1(n-1)=p1n-dx*(k11+2.d0*k12+2.d0*k13+k14)/6.d0
        q1(n-1)=q1n-dx*(m11+2.d0*m12+2.d0*m13+m14)/6.d0
        p2(n-1)=p2n-dx*(k21+2.d0*k22+2.d0*k23+k24)/6.d0
        q2(n-1)=q2n-dx*(m21+2.d0*m22+2.d0*m23+m24)/6.d0
        call dmag2op(pp1(n-1),qp1(n-1),pp2(n-1),qp2(n-1),
     >             p1(n-1),q1(n-1),p2(n-1),q2(n-1),
     >             xk1,xk2,facl1,facl2,cin,e,rnm1,
     >             vrnm1,ba11nm1,ba12nm1,ba22nm1,bb11nm1,bb22nm1)
c
      end do
c
c
c---> begin adams procedure (points ns+nout-6,  ... , 1)
c
      do n=ns+nout-5,2,-1
c
        x=xs+(n-ns)*dx
        r=dexp(x-dx)
        vrc=vr(n-1)
        ba11c=ba11(n-1)
        ba12c=ba12(n-1)
        ba22c=ba22(n-1)
        bb11c=bb11(n-1)
        bb22c=bb22(n-1)
c
        p1sum=646.d0*pp1(n)-264.d0*pp1(n+1)+106.d0*pp1(n+2)-
     >        19.d0*pp1(n+3)
        q1sum=646.d0*qp1(n)-264.d0*qp1(n+1)+106.d0*qp1(n+2)-
     >        19.d0*qp1(n+3)
        p2sum=646.d0*pp2(n)-264.d0*pp2(n+1)+106.d0*pp2(n+2)-
     >        19.d0*pp2(n+3)
        q2sum=646.d0*qp2(n)-264.d0*qp2(n+1)+106.d0*qp2(n+2)-
     >        19.d0*qp2(n+3)
c
c  predict for point n-1
c
        p10=p1(n)-dx*(251.d0*pp1(n+4)-1274.d0*pp1(n+3)+
     >      2616.d0*pp1(n+2)-2774.d0*pp1(n+1)+1901.d0*pp1(n))/720.d0
        q10=q1(n)-dx*(251.d0*qp1(n+4)-1274.d0*qp1(n+3)+
     >      2616.d0*qp1(n+2)-2774.d0*qp1(n+1)+1901.d0*qp1(n))/720.d0
        p20=p2(n)-dx*(251.d0*pp2(n+4)-1274.d0*pp2(n+3)+
     >      2616.d0*pp2(n+2)-2774.d0*pp2(n+1)+1901.d0*pp2(n))/720.d0
        q20=q2(n)-dx*(251.d0*qp2(n+4)-1274.d0*qp2(n+3)+
     >      2616.d0*qp2(n+2)-2774.d0*qp2(n+1)+1901.d0*qp2(n))/720.d0
c
c  correct
c
        if(imode.eq.1) then
c
        do nit=1,nitmax
c
          call dmag2op(pp1nm1,qp1nm1,pp2nm1,qp2nm1,
     >               p10,q10,p20,q20,
     >               xk1,xk2,facl1,facl2,cin,e,r,
     >               vrc,ba11c,ba12c,ba22c,bb11c,bb22c)
          p11=p1(n)-dx*(251.d0*pp1nm1+p1sum)/720.d0
          q11=q1(n)-dx*(251.d0*qp1nm1+q1sum)/720.d0
          p21=p2(n)-dx*(251.d0*pp2nm1+p2sum)/720.d0
          q21=q2(n)-dx*(251.d0*qp2nm1+q2sum)/720.d0
c
c  compare predictor with corrector
c
          if(test*cdabs(p11-p10).gt.cdabs(p10)) goto 15
          if(test*cdabs(q11-q10).gt.cdabs(q10)) goto 15
          if(test*cdabs(p21-p20).gt.cdabs(p20)) goto 15
          if(test*cdabs(q21-q20).gt.cdabs(q20)) goto 15
          goto 10
c
   15     p10=p11
          q10=q11
          p20=p21
          q20=q21
c
        end do
        write(6,*) n-1,r,nit,' not converged'
c
        else
c
        call dmag2op(pp1nm1,qp1nm1,pp2nm1,qp2nm1,
     >             p10,q10,p20,q20,
     >             xk1,xk2,facl1,facl2,cin,e,r,
     >             vrc,ba11c,ba12c,ba22c,bb11c,bb22c)
        p11=p1(n)-dx*(251.d0*pp1nm1+p1sum)/720.d0
        q11=q1(n)-dx*(251.d0*qp1nm1+q1sum)/720.d0
        p21=p2(n)-dx*(251.d0*pp2nm1+p2sum)/720.d0
        q21=q2(n)-dx*(251.d0*qp2nm1+q2sum)/720.d0
c
        q11=dkoef1*q11+dkoef2*q10
        p11=dkoef1*p11+dkoef2*p10
        q21=dkoef1*q21+dkoef2*q20
        p21=dkoef1*p21+dkoef2*p20
c
        end if
c
   10   q1(n-1)=q11
        p1(n-1)=p11
        q2(n-1)=q21
        p2(n-1)=p21
        call dmag2op(pp1(n-1),qp1(n-1),pp2(n-1),qp2(n-1),
     >             p1(n-1),q1(n-1),p2(n-1),q2(n-1),
     >             xk1,xk2,facl1,facl2,cin,e,r,
     >             vrc,ba11c,ba12c,ba22c,bb11c,bb22c)
c
      end do
c
c     store radial amplitudes times radius
c
      do n=1,ns+nout
        g(1,ii,n)=p1(n)
        g(2,ii,n)=p2(n)
        f(1,ii,n)=q1(n)/c
        f(2,ii,n)=q2(n)/c
      end do
      gp(1,ii)=(pp1(ns)-p1(ns))/rs
      gp(2,ii)=(pp2(ns)-p2(ns))/rs
c
      end do
c
      return
      end
      subroutine dmag2op(pp1,qp1,pp2,qp2,p1,q1,p2,q2,
     >                 xk1,xk2,facl1,facl2,cin,e,r,
     >                 vr,ba11,ba12,ba22,bb11,bb22)
c
      implicit real*8 (a-h,o-z)
      complex*16 pp1,qp1,pp2,qp2,p1,q1,p2,q2,e,tr,sr1,sr2
c
      tr=r*e-vr
      sr1=cin*tr+r+cin*bb11
      sr2=cin*tr+r+cin*bb22
c
      qp1= xk1*q1 - tr*p1 + ba11*p1 + ba12*p2 + facl1*p1/sr1
c
      pp1=-xk1*p1 + sr1*q1
c
      qp2= xk2*q2 - tr*p2 + ba22*p2 + ba12*p1 + facl2*p2/sr2
c
      pp2=-xk2*p2 + sr2*q2
c
      return
      end

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: hkershaw-brown/MuST


program first_program
!Este programa multiplica dos numeros enteros
implicit none
!declaracion de variables
integer :: A,B,C

write (*,*) 'Este programa multiplica dos numeros enteros'

print *, 'Escribe el valor de A'
read *, A

print *, 'Escribe el valor de B'
read *, B

C = A * B

print *, 'El producto de los dos numeros es:', C
print *, 'Man''s best friend'
print *, '"Who cares?"'
print *, 'Distance ='

end program


text:
module zone
!=======================================================================
! (formerly zone.h) global (3D) data arrays
!=======================================================================

 INTEGER, PARAMETER :: imax = 100, jmax = 100, kmax = 100    ! Memory dimensions
 INTEGER, PARAMETER :: pey = 4, pez = 4      ! number of MPI tasks
!              ####### for 2D:  ^^^  IF kmax=1, MUST HAVE pez=1   #############
 INTEGER, PARAMETER :: nvar = 6              ! number of primitive fluid variables

 INTEGER :: isy, isz, js, ks, Ya2abuff_size, Za2abuff_size
 INTEGER :: npe, npey, npez, mype, mypey, mypez  ! # of pes and local pe number
 INTEGER :: MPI_COMM_ROW, MPI_COMM_COL
 INTEGER :: VH1_DATATYPE

 INTEGER :: ngeomx, ngeomy, ngeomz       ! XYZ Geometry flag
 INTEGER :: nleftx, nlefty, nleftz       ! XYZ Lower Boundary Condition
 INTEGER :: nrightx,nrighty,nrightz      ! XYZ Upper Boundary Condition

 REAL, DIMENSION(imax,jmax/pey,kmax/pez) :: zro, zpr, zux, zuy, zuz, zfl
 REAL, DIMENSION(imax) :: zxa, zdx, zxc
 REAL, DIMENSION(jmax) :: zya, zdy, zyc
 REAL, DIMENSION(kmax) :: zza, zdz, zzc

 REAL, DIMENSION(nvar,kmax/pez,jmax/pey,imax) :: send1
 REAL, DIMENSION(nvar,kmax/pez,imax/pey,jmax) :: send2
 REAL, DIMENSION(nvar,jmax/pey,kmax/pez,imax) :: send3
 REAL, DIMENSION(nvar,jmax/pey,imax/pez,kmax) :: send4
 REAL, DIMENSION(nvar,kmax/pez,jmax/pey,imax/pey,pey) :: recv1
 REAL, DIMENSION(nvar,kmax/pez,imax/pey,jmax/pey,pey) :: recv2
 REAL, DIMENSION(nvar,jmax/pey,kmax/pez,imax/pez,pez) :: recv3
 REAL, DIMENSION(nvar,jmax/pey,imax/pez,kmax/pez,pez) :: recv4

 EQUIVALENCE ( send1(1,1,1,1),   send2(1,1,1,1),   send3(1,1,1,1),   send4(1,1,1,1) )
 EQUIVALENCE ( recv1(1,1,1,1,1), recv2(1,1,1,1,1), recv3(1,1,1,1,1), recv4(1,1,1,1,1) )

end module zone

max_stars_repo_name: samgeen/Weltgeist

max_issues_repo_name: samgeen/Weltgeist

max_forks_repo_name: samgeen/Weltgeist


text:
      REAL FUNCTION FMKIEEE(A)
C   $$$ SUBPROGRAM DOCUMENTATION BLOCK
C    . . . .
C    SUBPROGRAM: mkieee
C    PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-05-09
C
C    ABSTRACT: This subroutine stores a list of real values in
C    32-bit IEEE floating point format.
C
C    PROGRAM HISTORY LOG:
C    2000-05-09 Gilbert
C
C    USAGE: CALL mkieee(a)
C    INPUT ARGUMENT LIST:
C    a - Input floating point value.
C
C    OUTPUT ARGUMENT LIST: None.
C
C    REMARKS: None
C
C    ATTRIBUTES:
C    LANGUAGE: Fortran 90
C    MACHINE: IBM SP
C
C   $$$

      PARAMETER (two23=2.**23)
      PARAMETER (two126=2.**126)
C
      EQUIVALENCE(RTEMP,IEEE)
C
      alog2=alog(2.0)

      ieee=0

      if (a.ne.0.) then
C
C    Set Sign bit (bit 31 - leftmost bit)
C
         if (a.lt.0.0) then
           ieee=ibset(ieee,31)
           atemp=abs(a)
         else
           ieee=ibclr(ieee,31)
           atemp=a
         endif
C
C    Determine exponent n with base 2
C
         n=INT(flr(alog(atemp)/alog2))
         iexp=n+127
C    overflow
         if (n.gt.127) iexp=255
         if (n.lt.-127) iexp=0
C    set exponent bits ( bits 30-23 )
         call mvbits(iexp,0,8,ieee,23)
C
C    Determine Mantissa
C
         if (iexp.ne.255) then
           if (iexp.ne.0) then
             atemp=(atemp/(2.0**n))-1.0
           else
             atemp=atemp*two126
           endif
           imant=nint(atemp*two23)
         else
           imant=0
         endif
C    set mantissa bits ( bits 22-0 )
         call mvbits(imant,0,23,ieee,0)
      endif
C
      FMKIEEE=RTEMP
      return
      end

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
c
c     ****************************************************************
c     *                                                              *
c     *                 subroutine cohes_rot_mat                     *
c     *                                                              *
c     *                    written by : aroy                         *
c     *                 last modified : 02/10/13 rhd                 *
c     *                                                              *
c     *     compute rotation matrix bigR[span,3,3] that rotates      *
c     *     global vector quantities into the local orthogonal       *
c     *     system at the parametric center of interface elements.   *
c     *     the local system x-y axes lie in the tangent plane.      *
c     *     local z is normal to tangent plane                       *
c     *                                                              *
c     *     geometrically nonlinear elements: before entering this   *
c     *     subroutine, the ce vector has been modified such that    *
c     *     the first nnode/2 points describe the reference surface. *
c     *     see cohes_ref_surface.                                   *
c     *                                                              *
c     *                                                              *
c     ****************************************************************
c
      subroutine cohes_rot_mat(span, felem, nnode, etype, ce, bigR )
      implicit none
c
      include 'param_def'
c
c                  parameters
c
      integer :: span, felem, nnode, etype
      double precision :: ce(mxvl,mxecor), bigR(mxvl,3,3)
c
c                   local data
c
      integer :: i, j, imxvl, ispan, innode
      double precision ::
     & x21,y21,z21, x31,y31,z31, d1, d2, d3, a2, b2, c2,
     & a3, b3, c3, x1, y1, z1,
     & e123_to_ref(3,3),
     & ce_rotated(mxvl,mxecor), len, nzeta(mxndel),
     & neta(mxndel), nxi(mxndel), jac(mxvl,3,3), dj(mxvl),
     & g_to_e123(mxvl,3,3)
      double precision, parameter :: zero1 = 0.0d0, zero2 = 0.0d0,
     &                               zero3 = 0.0d0,
     &                               third1 = 1.0d0/3.0d0,
     &                               third2 = 1.0d0/3.0d0
c
      logical :: quad_interface, tri_interface
      logical, parameter :: local_debug = .false.
c
c
c               for geonl these are deformed coordinates of the
c               reference surface (most often the deformed
c               middle surface. see notes above and in routine
c               cohes_ref_surface.
c
c               Note: that first nnode/2 terms in
c               ce vector describe the reference surface.)
c
c               shift global element coords so element node 1
c               is at origin. really not necessary but makes
c               potential debugging easier.
c
      x1 = ce(1,1); y1 = ce(1,nnode+1); z1 = ce(1,2*nnode+1)
      do j = 1, nnode
!DIR$ IVDEP
         do i = 1, span
           ce(i,j) = ce(i,j) - x1
           ce(i,nnode+j) = ce(i,nnode+j) - y1
           ce(i,2*nnode+j) = ce(i,2*nnode+j) - z1
         end do
      end do
c
c               an element local system based on element nodes
c               1, 2,3. local z is normal to plane containing
c               element nodes 1, 2, 3. Then rotate the shifted
c               global coordinates into this 1-2-3 system.
c
!DIR$ IVDEP
      do i = 1, span
        x21 = ce(i,2) - ce(i,1)
        y21 = ce(i,nnode+2) - ce(i,nnode+1)
        z21 = ce(i,2*nnode+2) - ce(i,2*nnode+1)
        x31 = ce(i,3) - ce(i,1)
        y31 = ce(i,nnode+3) - ce(i,nnode+1)
        z31 = ce(i,2*nnode+3) - ce(i,2*nnode+1)
        d1  = sqrt(x21**2+y21**2+z21**2)
c
        g_to_e123(i,1,1) = x21/d1
        g_to_e123(i,1,2) = y21/d1
        g_to_e123(i,1,3) = z21/d1
c
        a3 = y21*z31 - y31*z21
        b3 = x31*z21 - x21*z31
        c3 = x21*y31 - x31*y21
        d3 = sqrt(a3**2+b3**2+c3**2)
c
        g_to_e123(i,3,1) = a3/d3
        g_to_e123(i,3,2) = b3/d3
        g_to_e123(i,3,3) = c3/d3
c
        a2 = g_to_e123(i,3,2)*g_to_e123(i,1,3) -
     &       g_to_e123(i,1,2)*g_to_e123(i,3,3)
        b2 = g_to_e123(i,1,1)*g_to_e123(i,3,3) -
     &       g_to_e123(i,3,1)*g_to_e123(i,1,3)
        c2 = g_to_e123(i,3,1)*g_to_e123(i,1,2) -
     &       g_to_e123(i,1,1)*g_to_e123(i,3,2)
        d2 = sqrt(a2**2+b2**2+c2**2)
c
        g_to_e123(i,2,1) = a2/d2
        g_to_e123(i,2,2) = b2/d2
        g_to_e123(i,2,3) = c2/d2
      end do
c
      imxvl = mxvl; ispan = span; innode = nnode
      call rotate_cohes_var( imxvl, ispan, innode, g_to_e123,
     &                       ce, ce_rotated )
c
c              we need the orthogonal system (s1, s2, n) at the
c              parametric center of the reference surface, where
c              s1 and s2 lie in the tangent plane to the reference
c              surface.
c
c              compute the orthogonal rotation R[3,3] that takes
c              vector values from the orthogonal 1-2-3 system
c              constructed above to s1, s2, n.
c
c              for quadrilateral interface elements, use (xi=eta=0)
c              for triangular interface elements, use s1=s2=1/3
c              (s3 is then also 1/3).
c
c              use rows 1,2 in computed [J] at parametric center
c              to define vectors (generally nonorthogonal) in tangent
c              plane. take cross product to get n-direction.
c              use additional cross-products to get orthogonal
c              vectors in tangent plane.
c
c              final coordinate transformation from global XYZ
c              system to interface system on and normal to tangent
c              plane [bigR] = [e123_to_ref] * [g_to_e123]
c
c              get derivates of shape functions at parametric center
c              and coordinate jacobian
c
       quad_interface = etype .eq. 12
       tri_interface  = etype .eq. 14  .or. etype .eq. 15
c
       nxi = zero1; neta = zero1; nzeta = zero1  ! prevents unint failure
       if( quad_interface )
     &   call derivs( etype, zero1, zero2, zero3, nxi, neta, nzeta )
c
       if( tri_interface )
     &   call derivs( etype, third1, third2, zero1, nxi, neta, nzeta )
c
       call jacob_cohes( etype, imxvl, ispan, innode, ce_rotated,
     &                   nxi, neta, nzeta, dj, jac, 2 )
c
!DIR$ IVDEP
       do i = 1, span
c
        d1 = sqrt(jac(i,1,1)*jac(i,1,1) + jac(i,1,2)*jac(i,1,2)
     &             +  jac(i,1,3)*jac(i,1,3) )
c
        e123_to_ref(1,1) = jac(i,1,1)/d1
        e123_to_ref(1,2) = jac(i,1,2)/d1
        e123_to_ref(1,3) = jac(i,1,3)/d1
c
        a3 = jac(i,1,2)*jac(i,2,3) - jac(i,2,2)*jac(i,1,3)
        b3 = jac(i,2,1)*jac(i,1,3) - jac(i,1,1)*jac(i,2,3)
        c3 = jac(i,1,1)*jac(i,2,2) - jac(i,2,1)*jac(i,1,2)
        d3 = sqrt(a3*a3+b3*b3+c3*c3)
c
        e123_to_ref(3,1) = a3/d3
        e123_to_ref(3,2) = b3/d3
        e123_to_ref(3,3) = c3/d3
c
        a2 = e123_to_ref(3,2)*e123_to_ref(1,3) -
     &       e123_to_ref(1,2)*e123_to_ref(3,3)
        b2 = e123_to_ref(1,1)*e123_to_ref(3,3) -
     &       e123_to_ref(3,1)*e123_to_ref(1,3)
        c2 = e123_to_ref(3,1)*e123_to_ref(1,2) -
     &       e123_to_ref(1,1)*e123_to_ref(3,2)
        d2 = sqrt(a2*a2+b2*b2+c2*c2)
c
        e123_to_ref(2,1) = a2/d2
        e123_to_ref(2,2) = b2/d2
        e123_to_ref(2,3) = c2/d2
c
        bigR(i,1,1) = e123_to_ref(1,1)*g_to_e123(i,1,1) +
     &                e123_to_ref(1,2)*g_to_e123(i,2,1) +
     &                e123_to_ref(1,3)*g_to_e123(i,3,1)
        bigR(i,1,2) = e123_to_ref(1,1)*g_to_e123(i,1,2) +
     &                e123_to_ref(1,2)*g_to_e123(i,2,2) +
     &                e123_to_ref(1,3)*g_to_e123(i,3,2)
        bigR(i,1,3) = e123_to_ref(1,1)*g_to_e123(i,1,3) +
     &                e123_to_ref(1,2)*g_to_e123(i,2,3) +
     &                e123_to_ref(1,3)*g_to_e123(i,3,3)
        bigR(i,2,1) = e123_to_ref(2,1)*g_to_e123(i,1,1) +
     &                e123_to_ref(2,2)*g_to_e123(i,2,1) +
     &                e123_to_ref(2,3)*g_to_e123(i,3,1)
        bigR(i,2,2) = e123_to_ref(2,1)*g_to_e123(i,1,2) +
     &                e123_to_ref(2,2)*g_to_e123(i,2,2) +
     &                e123_to_ref(2,3)*g_to_e123(i,3,2)
        bigR(i,2,3) = e123_to_ref(2,1)*g_to_e123(i,1,3) +
     &                e123_to_ref(2,2)*g_to_e123(i,2,3) +
     &                e123_to_ref(2,3)*g_to_e123(i,3,3)
        bigR(i,3,1) = e123_to_ref(3,1)*g_to_e123(i,1,1) +
     &                e123_to_ref(3,2)*g_to_e123(i,2,1) +
     &                e123_to_ref(3,3)*g_to_e123(i,3,1)
        bigR(i,3,2) = e123_to_ref(3,1)*g_to_e123(i,1,2) +
     &                e123_to_ref(3,2)*g_to_e123(i,2,2) +
     &                e123_to_ref(3,3)*g_to_e123(i,3,2)
        bigR(i,3,3) = e123_to_ref(3,1)*g_to_e123(i,1,3) +
     &                e123_to_ref(3,2)*g_to_e123(i,2,3) +
     &                e123_to_ref(3,3)*g_to_e123(i,3,3)
      end do
c
      if( local_debug ) then
        do i = 1, span
          write(*,100) bigR(i,1,1), bigR(i,1,2), bigR(i,1,3)
          write(*,100) bigR(i,2,1), bigR(i,2,2), bigR(i,2,3)
          write(*,100) bigR(i,3,1), bigR(i,3,2), bigR(i,3,3)
        end do
      end if
c
      return
c
100   format('Lambda ->',3(2X,f10.6))
101   format('R ->',3(2X,f10.6))
c
      end
c
c     ****************************************************************
c     *                                                              *
c     *             subroutine rotate_cohes_var                      *
c     *                                                              *
c     *                    written by : aroy                         *
c     *                    modified: 2/10/13 rhd                     *
c     *                                                              *
c     *           rotates multiple vectors stored in packed form     *
c     *                                                              *
c     ****************************************************************
c
      subroutine rotate_cohes_var( mxvl, span, nnode, rot, var, rvar )
      implicit none
c
c                  parameter declarations
c
      integer mxvl, span, nnode
      double precision
     & var(mxvl,*), rvar(mxvl,*), rot(mxvl,3,3)
c
c                  local declarations
c
      integer i, j, jj, k
c
c
c              var has nnode * 3 terms. 1-nnode are x components.
c              then the nnode y compoenents, then the nnode z
c              components.
c
      do i = 1, nnode
         j = i + nnode
        jj = i + 2*nnode
!DIR$ IVDEP
          do k = 1, span
             rvar(k,i) = var(k,i)*rot(k,1,1) + var(k,j)*rot(k,1,2)
     &                   + var(k,jj)*rot(k,1,3)
             rvar(k,j) = var(k,i)*rot(k,2,1) + var(k,j)*rot(k,2,2)
     &                   + var(k,jj)*rot(k,2,3)
             rvar(k,jj)= var(k,i)*rot(k,3,1) + var(k,j)*rot(k,3,2)
     &                   + var(k,jj)*rot(k,3,3)
         end do
      end do
c
      return
      end

c     ****************************************************************
c     *                                                              *
c     *                      subroutine jacob_cohes                  *
c     *                                                              *
c     *                       written by : aroy                      *
c     *                       last modified: 2/10/13 rhd             *
c     *                                                              *
c     *     compute 3x3 coordinate jacobian matrix & det at a point  *
c     *     for block of interface elements                          *
c     *                                                              *
c     ****************************************************************
c
c
c       uses the formula: [J] = [SFD] * [coords]
c         where   [J] = jacobian matrix
c               [SFD] = matrix of shape function derivatives
c                       (passed in as nxi, neta, nzeta)
c                [coords] = element nodal coordinates
c
      subroutine jacob_cohes( etype, mxvl, span, nnode, coords,
     &                        nxi, neta, nzeta, dj, jac, action )
      implicit none
c
c                  parameters
c
      integer etype, mxvl, span, nnode, action
c
      double precision
     & coords(mxvl,*), nzeta(*), neta(*), nxi(*), jac(mxvl,3,3), dj(*)
c
c                  locals
c
      integer i, j
      double precision
     & zero, j1, j2, j3
c
      data zero / 0.0d0 /
c
!DIR$ IVDEP
      do i = 1, span
         jac(i,1:3,1:3) = zero
         dj(i)          = zero
      end do
c
      do j = 1, nnode
!DIR$ IVDEP
         do i = 1, span
             jac(i,1,1)= jac(i,1,1)+nxi(j)*coords(i,j)
             jac(i,1,2)= jac(i,1,2)+nxi(j)*coords(i,nnode+j)
             jac(i,1,3)= jac(i,1,3)+nxi(j)*coords(i,2*nnode+j)
             jac(i,2,1)= jac(i,2,1)+neta(j)*coords(i,j)
             jac(i,2,2)= jac(i,2,2)+neta(j)*coords(i,nnode+j)
             jac(i,2,3)= jac(i,2,3)+neta(j)*coords(i,2*nnode+j)
             jac(i,3,1)= jac(i,3,1)+nzeta(j)*coords(i,j)
             jac(i,3,2)= jac(i,3,2)+nzeta(j)*coords(i,nnode+j)
             jac(i,3,3)= jac(i,3,3)+nzeta(j)*coords(i,2*nnode+j)
         end do
      end do
      if ( action .eq. 2 ) return
c
c           determinate of the jacobian matrix
c
!DIR$ IVDEP
      do i = 1, span
          j1= jac(i,2,2)*jac(i,3,3)-jac(i,2,3)*jac(i,3,2)
          j2= jac(i,2,1)*jac(i,3,3)-jac(i,2,3)*jac(i,3,1)
          j3= jac(i,2,1)*jac(i,3,2)-jac(i,2,2)*jac(i,3,1)
          dj(i)= jac(i,1,1)*j1-jac(i,1,2)*j2+jac(i,1,3)*j3
      end do
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine cohes_ref_surface                  *
c     *                                                              *
c     *                       written by : aroy                      *
c     *                                                              *
c     *                   last modified : 12/11/10 rhd comments      *
c     *                                                              *
c     *     for the geometrically nonlinear cohesive elements,       *
c     *     this subroutine computes the coordinates of the          *
c     *     reference surface for symmetric (top and bottom) and     *
c     *     non-specialized (middle) configurations.                 *
c     *     ce_upd is modified here to store the reference           *
c     *     surface coordinates.                                     *
c     *                                                              *
c     ****************************************************************
c
c          for cohesive elements, the strain-displacement matrix [B]
c          is given by [B] = [R][L][N]. for geometrically nonlinear
c          elements, the [R] matrix is computed using a
c          reference surface specified by the user. this
c          routine computes the reference surface coordinates
c          and stores them in ce_upd.
c
c          data structure in cu_upd:
c            on entry, cu_upd stores the current
c           (ie deformed) coords of element nodes in following form:
c
c   ce_upd = |        |
c            !   X    |\
c            | coor   |   \|    | node   1   of bottom surface      |
c            |        |   -'    | node   2   of bottom surface      |
c            |--------|         |    .............                  |
c            |        |         |    .............                  |
c            |   Y    |--->     | node 'nnode/2' of bottom surface  |
c            | coor   |         | node    1    of top    surface    |
c            |        |         | node    2    of top    surface    |
c            |--------|   -.    |    .............                  |
c            |        |   /|    |    .............                  |
c            |   Z    |  /      | node 'nnode/2' of top surface     |
c            | coor   | /
c            |        |
c
c           on exit from this routine, the first nnode/2 locations in
c           ce_upd store the coordinates of reference surface nodes.
c           the remaining nnode/2 locations store the reflected surface.
c           (the reflected surface is so called because at a later
c           point its coordinates will be reflected upon the
c           reference surface to take care of Poisson effect of bulk
c           material.)
c
c           the definition of reflected surface and the exact structure
c           of ce_upd on exit depends on reference surface (top, bottom
c           or middle). these are explained below in each case.
c
c           for symmetric models, the user must declare the "reference"
c           surface (top or bottom) as the surface that is connected t
c           o the solid elements. nodes on the reflected
c           surface will be given deformed coordinates
c           of the reference surface.
c
      subroutine cohes_ref_surface( span, mxvl, mxecor, surf, nnode,
     &                              totdof, ce, ce_upd, dj )
c
c              parameters
c
      implicit none
      integer span, mxvl, mxecor, surf, nnode, totdof
      double precision
     &    ce_upd(mxvl,mxecor), ce(mxvl,mxecor), dj(mxvl)
c
c              locals. one automatic array.
c
      integer i, j, k, shift, xb, xt, yb, yt, zb, zt
      double precision
     &  ce_refsurf(span,mxecor), zero, half
      data zero, half / 0.0d0, 0.5d0 /
      logical local_debug
      data local_debug /.false./
c
      shift = nnode/2
      xb = 0
      xt = xb + shift
      yb = nnode
      yt = yb + shift
      zb = 2*nnode
      zt = zb + shift
c
      if( local_debug ) then
        do i = 1, span
          write(*,*) 'surf = ',surf
        end do
      end if
c
      select case( surf )
c     ===================
c
      case( 1 )
c
c                 *** reference surface - top  ***
c
c      reference surface - top.
c      reflected surface - bottom.
c      change ce_upd such that on exit it has following structure.
c
c   ce_upd = |        |
c            |   X    |\
c            | coor   |   \|   | node   1   of reference (top) surface|
c            |        |   -'   | node   2   of reference (top) surface|
c            |--------|        |    .............                     |
c            |        |        |    .............                     |
c            |   Y    |---->   |node 'nnode/2' of reference (top) surf|
c            | coor   |        | node    1    of reflect (bot) surface|
c            |        |        | node    2    of reflect (bot) surface|
c            |--------|   -.   |    .............                     |
c            |        |   /|   |    .............                     |
c            |   Z    |  /     |node 'nnode/2' of reflect (bot) surf  |
c            | coor   | /
c            |        |
c
      do j = 1,shift
!DIR$ IVDEP
        do k = 1, span
            ce_refsurf(k,xb+j) = ce_upd(k,xt+j)
            ce_refsurf(k,xt+j) = ce_upd(k,xb+j)
            ce_refsurf(k,yb+j) = ce_upd(k,yt+j)
            ce_refsurf(k,yt+j) = ce_upd(k,yb+j)
            ce_refsurf(k,zb+j) = ce_upd(k,zt+j)
            ce_refsurf(k,zt+j) = ce_upd(k,zb+j)
         end do
      end do
c
      do j = 1, totdof
!DIR$ IVDEP
        ce_upd(1:span,j) = ce_refsurf(1:span,j)
      end do
c
      case( 2 )
c
c               *** reference surface - middle  ***
c
c      if penetration has occurred, define reference surface as
c      the original (undeformed) (top and bottom) surface of
c      cohesive element. the reflect surface is then defined by
c      0.5(deformed top+ deformed bottom).
c
c      if there is no penetration, define both reference and
c      reflect surface by 0.5(deformed top+ deformed bottom).
c
c
c      change ce_upd such that on exit it has following structure.
c
c   ce_upd = |        |
c            |  X     |\
c            | coor   |   \|   | node   1   of reference surface|
c            |        |   -'   | node   2   of reference surface|
c            |--------|        |    .............               |
c            |        |        |    .............               |
c            |   Y    |---->   |node 'nnode/2' of reference surf|
c            | coor   |        | node    1    of reflect surface|
c            |        |        | node    2    of reflect surface|
c            |--------|   -.   |    .............               |
c            |        |   /|   |    .............               |
c            |   Z    |  /     |node 'nnode/2' of reflect surf  |
c            | coor   | /
c            |        |
c
c
c
      do k = 1, span
      if( dj(k) .lt. zero ) then  ! interpenetration
             do j = 1, shift
               ce_refsurf(k,xb+j) =  ce(k,xb+j)
               ce_refsurf(k,xt+j) =
     &         half*( ce_upd(k,xb+j) + ce_upd(k,xt+j) )
               ce_refsurf(k,yb+j) =  ce(k,yb+j)
               ce_refsurf(k,yt+j) =
     &         half*( ce_upd(k,yb+j) + ce_upd(k,yt+j) )
               ce_refsurf(k,zb+j) =  ce(k,zb+j)
               ce_refsurf(k,zt+j) =
     &         half*( ce_upd(k,zb+j) + ce_upd(k,zt+j) )
            end do
      else
             do j = 1, shift
               ce_refsurf(k,xb+j) =
     &         half*( ce_upd(k,xb+j) + ce_upd(k,xt+j) )
               ce_refsurf(k,xt+j) =  ce_refsurf(k,xb+j)
               ce_refsurf(k,yb+j) =
     &         half*( ce_upd(k,yb+j) + ce_upd(k,yt+j) )
               ce_refsurf(k,yt+j) =  ce_refsurf(k,yb+j)
               ce_refsurf(k,zb+j) =
     &         half*( ce_upd(k,zb+j) + ce_upd(k,zt+j) )
               ce_refsurf(k,zt+j) =  ce_refsurf(k,zb+j)
            end do
        end if
      end do
c
      do j = 1, totdof
!DIR$ IVDEP
        ce_upd(1:span,j) = ce_refsurf(1:span,j)
      end do
c
      case( 3 )
c
c                *** reference surface - bottom  ***
c
c      reference surface - bottom.
c      reflect surface - top.
c      ce_upd already has following structure. no change required.
c
c   ce_upd = |        |
c            |   X    |\
c            | coor   |   \|   | node   1   of reference (bot) surface|
c            |        |   -'   | node   2   of reference (bot) surface|
c            |--------|        |    .............                     |
c            |        |        |    .............                     |
c            |   Y    |---->   |node 'nnode/2' of reference (bot) surf|
c            | coor   |        | node    1    of reflect (top) surface|
c            |        |        | node    2    of reflect (top) surface|
c            |--------|   -.   |    .............                     |
c            |        |   /|   |    .............                     |
c            |   Z    |  /     |node 'nnode/2' of reflect (top) surf  |
c            | coor   | /
c            |        |
c
c
      continue
c
      end select
      return
c
      end
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine cohes_mirror_refsurf               *
c     *                                                              *
c     *                  written by : aroy                           *
c     *                  last modified: 2/10/13  rhd                 *
c     *                                                              *
c     *     equates the updated nodal coordinates of the reference   *
c     *     surface to those of the reflected surface. this          *
c     *     is done only for a symmetric model and the reference     *
c     *     surface is either bottom or top. (In such cases, all     *
c     *     dof of the nodes attached to the reference plane must    *
c     *     be constrained)                                          *
c     *                                                              *
c     ****************************************************************
c
c           the reflected surface deforms with the bulk material.
c           thus its updated coordinates shows the Poisson's ratio
c           effect. however, the nodes attached to reference surface
c           have all dofs constrained. to obtain the true shape of the
c           deformed cohesive element the inplane coordinates of
c           reference surface should be equated to those of reflect
c           surface. here we equate all coordinates of these two
c           surfaces. since the interface cohesive element is
c           treated as a 2-D surface element (2x2 integration),
c           the equated normal coodinates will not introduce
c           any error. on the other hand, it makes the code
c           very simple :)
c
c            NOTE: before entering this subroutine, the ce vector has
c                been modified such that the first nnode/2 points
c                describe the reference surface and the remaining
c                nnode/2 points describe the reflect surface.
c                see SR cohes_ref_surface.
c
c            when the reference surface is middle this SR is a sham
c            one since the reference and reflect surface and their
c            updated coordinates are then the same by definition
c            (see  SR cohes_ref_surface).
c
      subroutine cohes_mirror_refsurf( span, mxvl, totdof,
     &    nnode, ce )
      implicit none
c
c                  parameters
c
      integer span, mxvl, totdof, nnode
      double precision
     & ce(mxvl,*)
c
c                  locals
c
      integer j, k, shift, xb, xt, yb, yt, zb, zt
c
      shift = nnode/2
c
      xb = 0
      xt = xb + shift
      yb = nnode
      yt = yb + shift
      zb = 2*nnode
      zt = zb + shift
c
      do k = 1, shift
       do j = 1, span
         ce(j,xb+k) = ce(j,xt+k)
         ce(j,yb+k) = ce(j,yt+k)
         ce(j,zb+k) = ce(j,zt+k)
       end do
      end do
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *             subroutine chk_cohes_penetrate                   *
c     *                                                              *
c     *                    written by : aroy                         *
c     *                last modified  : 02/10/13 rhd                 *
c     *                                                              *
c     *     support for subsequent check interpenetration of         *
c     *     interface-cohesive elements. checked using jacobian of   *
c     *     element with updated nodal coordinates. applicable only  *
c     *     to cohesive zone models which are isotropic under shear. *
c     *                                                              *
c     ****************************************************************
c
c           if penetration occurs, we need to set the
c           reference surface as the initial undeformed coordinates.
c           see SR cohes_ref_surf. coords passed in will be deformed
c           in current state.
c
      subroutine chk_cohes_penetrate( span, mxvl, felem, mxndel,
     &                                nnode, etype, coords, dj )
      implicit none
c
      integer span, mxvl, felem, mxndel, nnode, etype
      double precision
     &    coords(mxvl,*), dj(mxvl)
c
c              locally defined
c
      integer  solid_type
      logical  inter_quad_8, inter_tri_6, inter_tri_12
      double precision
     &    zero, half, third
c
c             automatic arrays
c
      double precision
     & nzeta(mxndel), neta(mxndel), nxi(mxndel), jac(mxvl,3,3)
       data zero, half, third / 0.0d0, 0.5d0, 0.333333333333333333d0 /
c
c         element interpenetration is checked later by the
c           sign (+ve or -ve) of the determinant of jacobian.
c
c         treat the interface elements as solids.
c           quadrilateral interface element -> brick
c           triangular interface element -> wedge
c
c        first compute shape function derivatives at element
c        parametric center.
c         for bricks, xi=eta=zeta=0.
c         for wedges, s1=s2=1/3, zeta=0.
c        (see shapef.f and derivs.f for triangle natural
c         coordinates s1, s2, s3.)
c
c        for quadrilateral cohesive element inter_8,
c        call derivs with etype=2 (8-node brick).
c        treat triangular elements as wedges.
c
      inter_quad_8 = etype .eq. 12
      inter_tri_6  = etype .eq. 14
      inter_tri_12 = etype .eq. 15
c
      if( inter_quad_8 ) then
         solid_type = 2
         call derivs( solid_type, zero, zero, zero, nxi, neta, nzeta )
      end if
c
      if(  inter_tri_6 )
     &    call deriv_trint6_wedge( etype, third, third, zero,
     &                             nxi, neta, nzeta)

      if( inter_tri_12 )
     &     call deriv_trint12_wedge( etype, third, third, zero,
     &                               nxi, neta, nzeta)
c
c           compute coordinate jacobian
c
       call jacob_cohes( etype, mxvl, span, nnode, coords,
     &                   nxi, neta, nzeta, dj, jac, 1 )
       return
       end
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine deriv_trint6_wedge                 *
c     *                                                              *
c     *                    written by : sushovan                     *
c     *                                 01/05/01                     *
c     *                                                              *
c     *     given a point in parametric coordinates, return the      *
c     *     parametric derivatives for each node shape function,     *
c     *     for the cohesive element trint6, treating it as wedge.   *
c     *                                                              *
c     ****************************************************************
c
c
c       use triangle natural coordinates (s1,s2,s3) for the wedge base.
c       (see SR shape14 and deriv14 for triangle natural coordinates).
c       zeta is along the height of the wedge. the two traingular faces
c       are at zeta = -1 (nodes 1,2,3) and zeta = +1 (nodes 4,5,6).
c
c       the shape functions for the wedge:
c         q1=0.5*s1*(1-zeta), q2=0.5*s2*(1-zeta), q3=0.5*s3*(1-zeta),
c         q4=0.5*s1*(1+zeta), q5=0.5*s2*(1+zeta), q6=0.5*s3*(1+zeta).
c
c       s3 is dependent coordinate. replace s3=1-s1-s2 while
c       taking the derivatives.
c
      subroutine deriv_trint6_wedge( etype, s1, s2, zeta, qs1, qs2,
     &                               qzeta )
      implicit none
c
      integer etype
      double precision
     &     s1, s2, zeta, qs1(*), qs2(*), qzeta(*)
c
c             local variables
c
      double precision
     &     s3, zero, half, one
      data zero, half, one /0.0d0, 0.5d0, 1.0d0/
c
      s3 = one - s1 -s2
c
c
c          dq_i/ds1, derivatives with respect to s1, the first
c          local coordinate : qs1(1) = dq1/ds1,
c                             qs1(2) = dq2/ds1,
c                             qs1(3) = dq3/ds1, etc.
c
      qs1(1) = half*(1-zeta)
      qs1(2) = zero
      qs1(3) = -half*(1-zeta)
      qs1(4) = half*(1+zeta)
      qs1(5) = zero
      qs1(6) = -half*(1+zeta)
c
c          dq_i/ds2, derivatives with respect to s2, the second
c          local coordinate: qs2(1) = dq1/ds2,
c                            qs2(2) = dq2/ds2,
c                            qs2(3) = dq3/ds2, etc.
c
      qs2(1) = zero
      qs2(2) = half*(1-zeta)
      qs2(3) = -half*(1-zeta)
      qs2(4) = zero
      qs2(5) = half*(1+zeta)
      qs2(6) = -half*(1+zeta)
c
c          dq_i/d(zeta), derivatives with respect to zeta, the third
c          local coordinate: qzeta(1) = dq1/d(zeta),
c                            qzeta(2) = dq2/d(zeta),
c                            qzeta(3) = dq3/d(zeta), etc.
c
      qzeta(1) = -half*s1
      qzeta(2) = -half*s2
      qzeta(3) = -half*s3
      qzeta(4) = half*s1
      qzeta(5) = half*s2
      qzeta(6) = half*s3
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine deriv_trint12_wedge                *
c     *                                                              *
c     *                    written by : sushovan                     *
c     *                                 01/05/01                     *
c     *                                                              *
c     *     given a point in parametric coordinates, return the      *
c     *     parametric derivatives for each node shape function,     *
c     *     for the cohesive element trint12, treating it as wedge.  *
c     *                                                              *
c     ****************************************************************
c
c
c       use triangle natural coordinates (s1,s2,s3) for the wedge base.
c       (see SR shape15 and deriv15 for triangle natural coordinates).
c       zeta is along the height of the wedge. the two traingular faces
c       are at zeta = -1 (nodes 1-6) and zeta = +1 (nodes 7-12).
c
c       shape functions for the wedge:
c
c        q1  = half*(2.0*s1*s1 - s1)(1-zeta)
c        q2  = half*(2.0*s2*s2 - s2)(1-zeta)
c        q3  = half*(2.0*s3*s3 - s3)(1-zeta)
c        q4  = half*four*s1*s2*(1-zeta)
c        q5  = half*four*s2*s3*(1-zeta)
c        q6  = half*four*s3*s1*(1-zeta)
c        q7  = half*(2.0*s1*s1 - s1)(1+zeta)
c        q8  = half*(2.0*s2*s2 - s2)(1+zeta)
c        q9  = half*(2.0*s3*s3 - s3)(1+zeta)
c        q10 = half*four*s1*s2*(1+zeta)
c        q11 = half*four*s2*s3*(1+zeta)
c        q12 = half*four*s3*s1*(1+zeta)
c
c       where half=0.5, four=4.0
c
c       s3 is dependent coordinate. replace s3=1-s1-s2 while
c       taking the derivatives.
c
      subroutine deriv_trint12_wedge( etype, s1, s2, zeta,
     &                                qs1, qs2, qzeta )
      implicit none
c
      integer etype
      double precision
     &     s1, s2, zeta, qs1(*), qs2(*), qzeta(*)
c
c             local variables
c
      double precision
     &     s3, zero, half, one, four
      data zero, half, one, four /0.0d0, 0.5d0, 1.0d0, 4.0d0/
c
      s3 = one - s1 -s2
c
c          dq_i/ds1, derivatives with respect to s1, the first
c          local coordinate : qs1(1) = dq1/ds1,
c                             qs1(2) = dq2/ds1,
c                             qs1(3) = dq3/ds1, etc.
c
      qs1(1)  = half*(four*s1 - one)*(1-zeta)
      qs1(2)  = zero
      qs1(3)  = half*(one - four*s3)*(1-zeta)
      qs1(4)  = half*four*s2*(1-zeta)
      qs1(5)  = -half*four*s2*(1-zeta)
      qs1(6)  = half*four*(s3 - s1)*(1-zeta)
      qs1(7)  = half*(four*s1 - one)*(1+zeta)
      qs1(8)  = zero
      qs1(9)  = half*(one - four*s3)*(1+zeta)
      qs1(10) = half*four*s2*(1+zeta)
      qs1(11) = -half*four*s2*(1+zeta)
      qs1(12) = half*four*(s3 - s1)*(1+zeta)
c
c          dq_i/ds2, derivatives with respect to s2, the second
c          local coordinate: qs2(1) = dq1/ds2,
c                            qs2(2) = dq2/ds2,
c                            qs2(3) = dq3/ds2, etc.
c
      qs2(1)  = zero
      qs2(2)  = half*(four*s2 - one)*(1-zeta)
      qs2(3)  = half*(one - four*s3)*(1-zeta)
      qs2(4)  = half*four*s1*(1-zeta)
      qs2(5)  = half*four*(s3 - s2)*(1-zeta)
      qs2(6)  = -half*four*s1*(1-zeta)
      qs2(7)  = zero
      qs2(8)  = half*(four*s2 - one)*(1+zeta)
      qs2(9)  = half*(one - four*s3)*(1+zeta)
      qs2(10) = half*four*s1*(1+zeta)
      qs2(11) = half*four*(s3 - s2)*(1+zeta)
      qs2(12) = -half*four*s1*(1+zeta)
c
c          dq_i/d(zeta), derivatives with respect to zeta, the third
c          local coordinate: qzeta(1) = dq1/d(zeta),
c                            qzeta(2) = dq2/d(zeta),
c                            qzeta(3) = dq3/d(zeta), etc.
c
      qzeta(1)  = -half*(2.0*s1*s1 - s1)
      qzeta(2)  = -half*(2.0*s2*s2 - s2)
      qzeta(3)  = -half*(2.0*s3*s3 - s3)
      qzeta(4)  = -half*four*s1*s2
      qzeta(5)  = -half*four*s2*s3
      qzeta(6)  = -half*four*s3*s1
      qzeta(7)  = half*(2.0*s1*s1 - s1)
      qzeta(8)  = half*(2.0*s2*s2 - s2)
      qzeta(9)  = half*(2.0*s3*s3 - s3)
      qzeta(10) = half*four*s1*s2
      qzeta(11) = half*four*s2*s3
      qzeta(12) = half*four*s3*s1
c
      return
      end

c
c
c     ****************************************************************
c     *                                                              *
c     *                 subroutine check_cohes_tri                   *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                   last modified: 2/18/13 rhd                 *
c     *                                                              *
c     *     make some sanity checks on the geometry of triangular    *
c     *     interface elements. uses undeformed coordinates.         *
c     *     prints warning messages if potential issues are found    *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine check_cohes_tri( etype, span, mxvl, felem, nnode,
     &                            iout,  ce )
      implicit none
c
      integer etype, span, mxvl, felem, nnode, iout
      double precision
     & ce(mxvl,*)
c
c                 local data. some are automatic arrays.
c
      double precision
     & xb(nnode), yb(nnode), zb(nnode), xt(nnode), yt(nnode),
     & zt(nnode), dxbt, dybt, dzbt, lbt,
     & x(nnode), y(nnode), z(nnode), xm, ym, zm, dx, dy, dz, dmid,
     & dx12, dy12, dz12, dx23, dy23, dz23, dx13, dy13, dz13,
     & lb12, lb23, lb13, lt12, lt23, lt13, l13, tolmin, tol_edge,
     & long_edge(span), short_edge, mid_edge_tol, top_bott_tol,
     & half
c
      logical flg_bottom, flg_top, collapsed

      integer nsurf, i, j, trint12_list(3,6), n1, n2, n3
      data tolmin, tol_edge, mid_edge_tol, top_bott_tol, half
     &  / 1.0d-08, 1.0d-03, 0.2d0, 0.05d0, 0.5d0 /
      data trint12_list / 1,4,2,  2,5,3,  3,6,1,  7,10,8,  8,11,9,
     &                    9,12,7 /
c
c             1. all 3 sides must have lengths > 0.0 (tolmin)
c             2. no side can have less length tol_edge * longest side
c
c             These checks treat both triangle interface elements as
c             six noded.
c
      nsurf = nnode / 2
      collapsed = .false.
      do i = 1, span
!DIR$ IVDEP
         do j = 1, nsurf
           xb(j) = ce(i,j)
           yb(j) = ce(i,nnode+j)
           zb(j) = ce(i,2*nnode+j)
           xt(j) = ce(i,nsurf+j)
           yt(j) = ce(i,nnode+nsurf+j)
           zt(j) = ce(i,2*nnode+nsurf+j)
         end do
         dx12 = xb(2) - xb(1)
         dy12 = yb(2) - yb(1)
         dz12 = zb(2) - zb(1)
         dx23 = xb(3) - xb(2)
         dy23 = yb(3) - yb(2)
         dz23 = zb(3) - zb(2)
         dx13 = xb(3) - xb(1)
         dy13 = yb(3) - yb(1)
         dz13 = zb(3) - zb(1)
         lb12 = sqrt( dx12*dx12 + dy12*dy12 + dz12*dz12 )
         lb23 = sqrt( dx23*dx23 + dy23*dy23 + dz23*dz23 )
         lb13 = sqrt( dx13*dx13 + dy13*dy13 + dz13*dz13 )
         flg_bottom = lb12 .le. tolmin .or. lb23 .le. tolmin .or.
     &       lb13 .le. tolmin
         dx12 = xt(2) - xt(1)
         dy12 = yt(2) - yt(1)
         dz12 = zt(2) - zt(1)
         dx23 = xt(3) - xt(2)
         dy23 = yt(3) - yt(2)
         dz23 = zt(3) - zt(2)
         dx13 = xt(3) - xt(1)
         dy13 = yt(3) - yt(1)
         dz13 = zt(3) - zt(1)
         lt12 = sqrt( dx12*dx12 + dy12*dy12 + dz12*dz12 )
         lt23 = sqrt( dx23*dx23 + dy23*dy23 + dz23*dz23 )
         lt13 = sqrt( dx13*dx13 + dy13*dy13 + dz13*dz13 )
         flg_top = lt12 .le. tolmin .or. lt23 .le. tolmin .or.
     &       lt13 .le. tolmin
         long_edge(i) = max( lb12, lb23, lb13, lt12, lt23, lt13 )
         short_edge = min( lb12, lb23, lb13, lt12, lt23, lt13 )
         if( short_edge .le. tol_edge * long_edge(i) )
     &         collapsed = .true.
         if( flg_bottom .or. flg_top .or. collapsed ) then
            write(iout,9000) felem+i-1
         end if
      end do
c
c             trint12 elements. check that location of mid-side
c             nodes are not too far from mid-point on straight
c             line connecting two corners.

      if( nnode .ne. 12 ) go to 300
      do i = 1, span
!DIR$ IVDEP
         do j = 1, 12
           x(j) = ce(i,j)
           y(j) = ce(i,12+j)
           z(j) = ce(i,24+j)
         end do
!DIR$ IVDEP
         do j = 1, 6
           n1 = trint12_list(1,j)
           n2 = trint12_list(2,j)  ! middle node on edge
           n3 = trint12_list(3,j)
           dx13 = x(n3) - x(n1)
           dy13 = y(n3) - y(n1)
           dz13 = z(n3) - z(n1)
           l13 = sqrt( dx13*dx13 + dy13*dy13 + dz13*dz13 )
           xm = x(n1) + half * dx13
           ym = y(n1) + half * dy13
           zm = z(n1) + half * dz13
           dx = x(n2) - xm
           dy = y(n2) - ym
           dz = z(n2) - zm
           dmid = sqrt( dx*dx + dy*dy + dz*dz )
           if( dmid .gt. mid_edge_tol * l13 ) then
             write(iout,9100) felem+i-1, n2
           end if
         end do
      end do
c
c             check distance between each pair of top and bottom
c             nodes

 300  continue
      nsurf = nnode / 2
      do i = 1, span
!DIR$ IVDEP
         do j = 1, nsurf
           xb(j) = ce(i,j)
           yb(j) = ce(i,nnode+j)
           zb(j) = ce(i,2*nnode+j)
           xt(j) = ce(i,nsurf+j)
           yt(j) = ce(i,nnode+nsurf+j)
           zt(j) = ce(i,2*nnode+nsurf+j)
           dxbt = xt(j) - xb(j)
           dybt = yt(j) - yb(j)
           dzbt = zt(j) - zb(j)
           lbt  = sqrt( dxbt*dxbt + dybt*dybt + dzbt*dzbt )
           if( lbt .gt. top_bott_tol * long_edge(i) ) then
            write(iout,9200) felem+i-1, j, j+nsurf, lbt, long_edge(i)
           end if
         end do
      end do

       return
c
 9000 format(/1x,'>>>>> Warning: interface-cohesive element: ',i8,
     & /,1x,'               appears to have one or more collapsed',
     &      ' edges.',//)
 9100 format(/1x,'>>>>> Warning: interface-cohesive element: ',i8,
     & /,1x,    '      Element mid-side node: ',i3,' appears to be',
     &          ' located a significant distance',
     & /,1x,    '      from the mid-point on a line connecting the',
     &          ' two corner nodes.',
     & /,1x,    '      Check coordinates & incidences for errors.',//)
 9200 format(/1x,'>>>>> Warning: interface-cohesive element: ',i8,
     & /,1x,    '      Element bottom-top pair node: ',2i3,' are',
     & /,1x,    '      are a significant distance apart: ',e14.6,
     & /,1x,    '      compared to longest edge length of element: ',
     & e14.6,//)
       end
c
c
c     ****************************************************************
c     *                                                              *
c     *                 subroutine check_cohes_quad                  *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                   last modified: 2/18/13 rhd                 *
c     *                                                              *
c     *     make some sanity checks on the geometry of the 8-node    *
c     *     interface element. uses undeformed coordinates.          *
c     *     prints warning messages if potential issues are found    *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine check_cohes_quad( etype, span, mxvl, felem, nnode,
     &                             iout, ce )
      implicit none
c
      integer etype, span, mxvl, felem, nnode, iout
      double precision
     & ce(mxvl,*)
c
c                 local data. some are automatic arrays.
c
      double precision
     & xb(nnode), yb(nnode), zb(nnode), xt(nnode), yt(nnode),
     & zt(nnode), dxbt, dybt, dzbt, lbt,
     & x(nnode), y(nnode), z(nnode), dx41, dy41, dz41,
     & dx12, dy12, dz12, dx23, dy23, dz23,
     & dx34, dy34, dz34, lt34, lb34,
     & lb12, lb23, lb41, lt12, lt23, lt41, tolmin, tol_edge,
     & long_edge(span), short_edge, top_bott_tol
c
      logical flg_bottom, flg_top, collapsed

      integer nsurf, i, j
      data tolmin, tol_edge, top_bott_tol
     &  / 1.0d-08, 1.0d-03, 0.05d0 /
c
c             1. all 4 sides must have lengths > 0.0 (tolmin)
c             2. no side can have less length tol_edge * longest side
c
      nsurf = nnode / 2
      collapsed = .false.
      do i = 1, span
!DIR$ IVDEP
         do j = 1, nsurf
           xb(j) = ce(i,j)
           yb(j) = ce(i,nnode+j)
           zb(j) = ce(i,2*nnode+j)
           xt(j) = ce(i,nsurf+j)
           yt(j) = ce(i,nnode+nsurf+j)
           zt(j) = ce(i,2*nnode+nsurf+j)
         end do
c
         dx12 = xb(2) - xb(1)
         dy12 = yb(2) - yb(1)
         dz12 = zb(2) - zb(1)
c
         dx23 = xb(3) - xb(2)
         dy23 = yb(3) - yb(2)
         dz23 = zb(3) - zb(2)
c
         dx34 = xb(4) - xb(3)
         dy34 = yb(4) - yb(3)
         dz34 = zb(4) - zb(3)
c
         dx41 = xb(1) - xb(4)
         dy41 = yb(1) - yb(4)
         dz41 = zb(1) - zb(4)
c
         lb12 = sqrt( dx12*dx12 + dy12*dy12 + dz12*dz12 )
         lb23 = sqrt( dx23*dx23 + dy23*dy23 + dz23*dz23 )
         lb34 = sqrt( dx34*dx34 + dy34*dy34 + dz34*dz34 )
         lb41 = sqrt( dx41*dx41 + dy41*dy41 + dz41*dz41 )
         flg_bottom = lb12 .le. tolmin .or. lb23 .le. tolmin .or.
     &       lb34 .le. tolmin  .or. lb41 .le. tolmin
c
         dx12 = xt(2) - xt(1)
         dy12 = yt(2) - yt(1)
         dz12 = zt(2) - zt(1)
c
         dx23 = xt(3) - xt(2)
         dy23 = yt(3) - yt(2)
         dz23 = zt(3) - zt(2)
c
         dx34 = xt(4) - xt(3)
         dy34 = yt(4) - yt(3)
         dz34 = zt(4) - zt(3)
c
         dx41 = xt(1) - xt(4)
         dy41 = yt(1) - yt(4)
         dz41 = zt(1) - zt(4)
c
         lt12 = sqrt( dx12*dx12 + dy12*dy12 + dz12*dz12 )
         lt23 = sqrt( dx23*dx23 + dy23*dy23 + dz23*dz23 )
         lt34 = sqrt( dx34*dx34 + dy34*dy34 + dz34*dz34 )
         lt41 = sqrt( dx41*dx41 + dy41*dy41 + dz41*dz41 )
         flg_top = lt12 .le. tolmin .or. lt23 .le. tolmin .or.
     &          lt34 .le. tolmin  .or. lt41 .le. tolmin
c
         long_edge(i) = max( lb12, lb23, lb34, lb41, lt12, lt23,
     &                       lt34, lt41)
         short_edge = min( lb12, lb23, lb34, lb41, lt12, lt23,
     &                       lt34, lt41 )
         if( short_edge .le. tol_edge * long_edge(i) )
     &         collapsed = .true.
         if( flg_bottom .or. flg_top .or. collapsed ) then
            write(iout,9000) felem+i-1
         end if
      end do
c
c             check distance between each pair of top and bottom
c             nodes
c
      nsurf = nnode / 2
      do i = 1, span
!DIR$ IVDEP
         do j = 1, nsurf
           xb(j) = ce(i,j)
           yb(j) = ce(i,nnode+j)
           zb(j) = ce(i,2*nnode+j)
           xt(j) = ce(i,nsurf+j)
           yt(j) = ce(i,nnode+nsurf+j)
           zt(j) = ce(i,2*nnode+nsurf+j)
           dxbt = xt(j) - xb(j)
           dybt = yt(j) - yb(j)
           dzbt = zt(j) - zb(j)
           lbt  = sqrt( dxbt*dxbt + dybt*dybt + dzbt*dzbt )
           if( lbt .gt. top_bott_tol * long_edge(i) ) then
            write(iout,9200) felem+i-1, j, j+nsurf, lbt, long_edge(i)
           end if
         end do
      end do

       return
c
 9000 format(/1x,'>>>>> Warning: interface-cohesive element: ',i8,
     & /,1x,'               appears to have one or more collapsed',
     &      ' edges.',//)
 9200 format(/1x,'>>>>> Warning: interface-cohesive element: ',i8,
     & /,1x,    '      Element bottom-top pair node: ',2i3,' are',
     & /,1x,    '      are a significant distance apart: ',e14.6,
     & /,1x,    '      compared to longest edge length of element: ',
     & e14.6,//)
c
       end