SUBROUTINE stiff4(km,coord,ym,pr)
!
! This subroutine generates the "analytical" stiffness matrix
! for a four node quadrilateral in plane strain based on nip=4.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::coord(:,:),ym,pr
 REAL(iwp),INTENT(OUT)::km(:,:)
 REAL(iwp)::e1,e2,g,x1,x2,x3,x4,y1,y2,y3,y4,a2,a2st3,alph,beta,c11,c21,   &
   c31,c41,c51,c61,c71,c81,c22,c32,c42,c52,c62,c72,c82,c33,c43,c53,c63,   &
   c73,c83,c44,c54,c64,c74,c84,c55,c65,c75,c85,c66,c76,c86,c77,c87,c88,   &
   f1,f2,s1,s2,s3,s4,t1,t2,t3,t4,pt5=0.5_iwp,one=1.0_iwp,two=2.0_iwp,     &
   d3=3.0_iwp,zero=0.0_iwp
 INTEGER::i,j
 e1=ym*(one-pr)/(one+pr)/(1-two*pr)
 e2=pr*e1/(one-pr)
 g=ym/two/(one+pr)
!
 x1=coord(1,1)
 x2=coord(2,1)
 x3=coord(3,1)
 x4=coord(4,1)
 y1=coord(1,2)
 y2=coord(2,2)
 y3=coord(3,2)
 y4=coord(4,2)
!
 a2=(x4-x2)*(y3-y1)-(x3-x1)*(y4-y2)
 a2st3=d3*a2*a2
!
 CALL groupa(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c11=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(1,1)=c11
!
 CALL groupa(x2,x3,x4,x1,y2,y3,y4,y1,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c33=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(3,3)=c33
!
 CALL groupa(x3,x4,x1,x2,y3,y4,y1,y2,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c55=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(5,5)=c55
!
 CALL groupa(x4,x1,x2,x3,y4,y1,y2,y3,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c77=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(7,7)=c77
!
 CALL groupa(y4,y3,y2,y1,x4,x3,x2,x1,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c88=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,8)=c88
!
 CALL groupa(y1,y4,y3,y2,x1,x4,x3,x2,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c22=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(2,2)=c22
!
 CALL groupa(y2,y1,y4,y3,x2,x1,x4,x3,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c44=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(4,4)=c44
!
 CALL groupa(y3,y2,y1,y4,x3,x2,x1,x4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c66=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(6,6)=c66
!
 CALL groupb(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c21=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(2,1)=c21
!
 CALL groupb(x2,x3,x4,x1,y2,y3,y4,y1,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c43=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(4,3)=c43
!
 CALL groupb(x3,x4,x1,x2,y3,y4,y1,y2,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c65=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(6,5)=c65
!
 CALL groupb(x4,x1,x2,x3,y4,y1,y2,y3,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c87=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,7)=c87
!
 CALL groupc(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c31=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(3,1)=c31
!
 CALL groupc(x2,x3,x4,x1,y2,y3,y4,y1,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c53=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(5,3)=c53
!
 CALL groupc(x3,x4,x1,x2,y3,y4,y1,y2,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c75=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(7,5)=c75
!
 CALL groupc(x4,x1,x2,x3,y4,y1,y2,y3,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c71=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(7,1)=c71
!
 CALL groupc(y4,y3,y2,y1,x4,x3,x2,x1,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c86=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,6)=c86
!
 CALL groupc(y1,y4,y3,y2,x1,x4,x3,x2,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c82=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,2)=c82
!
 CALL groupc(y2,y1,y4,y3,x2,x1,x4,x3,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c42=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(4,2)=c42
!
 CALL groupc(y3,y2,y1,y4,x3,x2,x1,x4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c64=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(6,4)=c64
!
 CALL groupd(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c41=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(4,1)=c41
!
 CALL groupd(x2,x3,x4,x1,y2,y3,y4,y1,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c63=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(6,3)=c63
!
 CALL groupd(x3,x4,x1,x2,y3,y4,y1,y2,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x3,x4,x1,x2,y3,y4,y1,y2,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c85=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,5)=c85
!
 CALL groupd(x4,x1,x2,x3,y4,y1,y2,y3,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x4,x1,x2,x3,y4,y1,y2,y3,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c72=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(7,2)=c72
!
 CALL groupd(y1,y2,y3,y4,x1,x2,x3,x4,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c32=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(3,2)=c32
!
 CALL groupd(y2,y3,y4,y1,x2,x3,x4,x1,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c54=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(5,4)=c54
!
 CALL groupd(y3,y4,y1,y2,x3,x4,x1,x2,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x3,x4,x1,x2,y3,y4,y1,y2,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c76=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(7,6)=c76
!
 CALL groupd(y4,y1,y2,y3,x4,x1,x2,x3,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x4,x1,x2,x3,y4,y1,y2,y3,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c81=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,1)=c81
!
 CALL groupe(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c51=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(5,1)=c51
!
 CALL groupe(x2,x3,x4,x1,y2,y3,y4,y1,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c73=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(7,3)=c73
!
 CALL groupe(y2,y1,y4,y3,x2,x1,x4,x3,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c84=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,4)=c84
!
 CALL groupe(y3,y2,y1,y4,x3,x2,x1,x4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c62=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(6,2)=c62
!
 CALL groupf(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c61=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(6,1)=c61
!
 CALL groupf(x2,x3,x4,x1,y2,y3,y4,y1,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c83=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(8,3)=c83
!
 CALL groupf(y1,y2,y3,y4,x1,x2,x3,x4,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c52=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(5,2)=c52
!
 CALL groupf(y2,y3,y4,y1,x2,x3,x4,x1,s1,s2,s3,s4,t1,t2,t3,t4)
 CALL f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c74=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*pt5
 km(7,4)=c74
!
 DO i=1,8
   DO j=i+1,8
     km(i,j)=km(j,i)
   END DO
 END DO
!
RETURN
CONTAINS
SUBROUTINE groupa(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::s1,s2,s3,s4,t1,t2,t3,t4,f1,f2
 s1=two*(y4-y2)**2
 s2=two*(x4-x2)**2
 s3=-s1/two
 s4=-s2/two
 t1=(y2-y3)**2+(y3-y4)**2+(y4-y2)**2
 t2=(x2-x3)**2+(x3-x4)**2+(x4-x2)**2
 t3=(y4-y3)**2-(y3-y2)**2
 t4=(x4-x3)**2-(x3-x2)**2
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-two*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-two*(x3*y1-x1*y3)
RETURN
END SUBROUTINE groupa
SUBROUTINE groupb(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::s1,s2,s3,s4,t1,t2,t3,t4,f1,f2
 s1=two*(x2-x4)*(y4-y2)
 s2=s1
 s3=-s1/2
 s4=s3
 t1=x2*(y4-two*y2+y3)+x3*(y2-two*y3+y4)+x4*(y2-two*y4+y3)
 t2=t1
 t3=x2*(y2-y3)+x3*(y4-y2)+x4*(y3-y4)
 t4=t3
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-two*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-two*(x3*y1-x1*y3)
RETURN
END SUBROUTINE groupb
SUBROUTINE groupc(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::s1,s2,s3,s4,t1,t2,t3,t4,f1,f2
 s1=(y4-y2)*(two*y1-y3-y4)
 s2=(x4-x2)*(two*x1-x3-x4)
 s3=(y4-y2)*(y4-y1)
 s4=(x4-x2)*(x4-x1)
 t1=(y3-y1)*(two*y2-y3-y4)
 t2=(x3-x1)*(two*x2-x3-x4)
 t3=(y3-y1)*(y3-y2)
 t4=(x3-x1)*(x3-x2)
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-two*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-two*(x3*y1-x1*y3)
RETURN
END SUBROUTINE groupc
SUBROUTINE groupd(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::s1,s2,s3,s4,t1,t2,t3,t4
 s1=(x3-x1)*(y4-y2)+(x4-x1)*(y4-y2)
 s2=(y3-y1)*(x4-x2)+(y4-y1)*(x4-x2)
 s3=(x4-x1)*(y2-y4)
 s4=(y4-y1)*(x2-x4)
 t1=(x3-x1)*(y4-y2)+(x3-x1)*(y3-y2)
 t2=(y3-y1)*(x4-x2)+(y3-y1)*(x3-x2)
 t3=(x3-x1)*(y2-y3)
 t4=(y3-y1)*(x2-x3)
RETURN
END SUBROUTINE groupd
SUBROUTINE groupe(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::s1,s2,s3,s4,t1,t2,t3,t4,f1,f2
 s1=-(y4-y2)**2
 s2=-(x4-x2)**2
 s3=zero
 s4=zero
 t1=(y3+y1)*(y4+y2)-two*(y4-y2)**2-two*(y1*y3+y2*y4)
 t2=(x3+x1)*(x4+x2)-two*(x4-x2)**2-two*(x1*x3+x2*x4)
 t3=(y4-y2)*(y1-y2+y3-y4)
 t4=(x4-x2)*(x1-x2+x3-x4)
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-two*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-two*(x3*y1-x1*y3)
RETURN
END SUBROUTINE groupe
SUBROUTINE groupf(x1,x2,x3,x4,y1,y2,y3,y4,s1,s2,s3,s4,t1,t2,t3,t4)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::s1,s2,s3,s4,t1,t2,t3,t4
 s1=(x4-x2)*(y4-y2)
 s2=s1
 s3=zero
 s4=zero
 t1=(x4-x2)*(y4-y2)+(x2-x1)*(y2-y3)+(x4-x1)*(y4-y3)
 t2=(y4-y2)*(x4-x2)+(y2-y1)*(x2-x3)+(y4-y1)*(x4-x3)
 t3=(x2-x1)*(y3-y2)+(x4-x1)*(y4-y3)
 t4=(y2-y1)*(x3-x2)+(y4-y1)*(x4-x3)
RETURN
END SUBROUTINE groupf
SUBROUTINE f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::f1,f2
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-two*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-two*(x3*y1-x1*y3)
RETURN
END SUBROUTINE f1f2
END SUBROUTINE stiff4


module derivatives
  use nrtype
contains
  subroutine constant_accel(x,y,dydx)
    implicit none
    real(dp), intent(in)  :: x, y(:)
    real(dp), intent(out) :: dydx(:)
    dydx = -9.81
  end subroutine constant_accel
end module derivatives


text:
great_ne(dd1,y1).
great_ne(hh1,b1).
great_ne(z1,aa1).
great_ne(u1,w1).
great_ne(u1,ff1).
great_ne(dd1,w1).
great_ne(kk1,h1).
great_ne(p1,a1).
great_ne(z1,y1).
great_ne(p1,c1).
great_ne(ii1,y1).
great_ne(u1,ee1).
great_ne(q1,ii1).
great_ne(q1,v1).
great_ne(z1,h1).
great_ne(p1,dd1).
great_ne(jj1,v1).
great_ne(q1,w1).
great_ne(kk1,i1).
great_ne(x1,b1).
great_ne(p1,b1).
great_ne(k1,ii1).
great_ne(cc1,v1).
great_ne(q1,ee1).
great_ne(z1,hh1).
great_ne(u1,dd1).
great_ne(t1,w1).
great_ne(f1,c1).
great_ne(t1,jj1).
great_ne(cc1,u1).
great_ne(bb1,h1).
great_ne(q1,a1).
great_ne(z1,t1).
great_ne(jj1,bb1).
great_ne(ee1,i1).
great_ne(q1,ff1).
great_ne(ee1,y1).

ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
c========================================================================
c
c     subroutine OPVAR
c
c     Computing optimization for ALLOCVAR - Value-at-Risk constraint
c
c-----------------------------------------------------------------------
      SUBROUTINE opvar ( n, cov, rho, neq, nin, ccst, bcst,
     &                   cinf, csup, vrimax, conflp, maxdic, epsdic,
     &                   iwork, dwork, wopt, vriopt, info )
c-----------------------------------------------------------------------
c
c     INPUT
c            n      : number of assets                           integer
c            cov    : covariance matrix (n*n)                     double
c            rho    : mean returns vector (n)                     double
c            neq    : number of initial equality constraints     integer
c            nin    : number of initial inequality constraints   integer
c            ccst   : matrix of constraints (n*(neq+nin))         double
c            bcst   : vector of constraints  vector(neq+nin)      double
c            cinf   : lower bounds (n)                            double
c            csup   : upper bounds (n)                            double
c            vrimax : Value-at-Risk constraint                    double
c            conflp : confidence level parameter ( 0<conflp<1 )   double
c            maxdic : maximum of dichotomy iterations            integer
c            epsdic : precision of dichotomy                      double
c
c     WORKSPACE
c            iwork  : 3*n + neq + 2*nin + 3                      integer
c            dwork  : n*(n + neq + nin + 10) + 2*neq + 4*nin + 4  double
c
c     OUTPUT
c            wopt   : optimal portfolio (n)                       double
c            vriopt : optimal Value-at-Risk                       double
c            info   : diagnostic argument                        integer
c
c     CALL
c            EVBOR  : computing the minimum and maximum of a vector
c            OPMV   : computing optimization for ALLOCMV
c            WZERO  : put zero at the values of a vector less than eps
c                     and increase the greatest element with residual
c            YV     : copy a vector in a vector
c            XV     : computing scalar product of two vectors = V1'*V2
c            OVTMCV : computing V'*M*V = scalar
c                    ( M vectorized square matrix(n*n), V vector(n) )
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, neq, nin, maxdic, info
      DOUBLE PRECISION vrimax, conflp, epsdic, vriopt
      DOUBLE PRECISION cov(*), rho(*), cinf(*), csup(*),
     &                 ccst(*), bcst(*), wopt(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER niter, infmax, infqua, piwo, pdwopt, pdwo
      DOUBLE PRECISION p, q, zc, rhopt, optmean, optvar, volat, vrisav,
     &                 rhomin, rhomax, epsvri, xvri, eps, ornumber
      PARAMETER ( eps = 1.0E-15 )
c
c     external subroutines
      EXTERNAL EVBOR, opmv, WZERO, YV, XV, OVTMCV, dinvnr
      DOUBLE PRECISION dinvnr
c
c     intrinsic functions
      INTRINSIC sqrt, min, max, abs
c
c-----------------------------------------------------------------------
c
c     initializations
      info = 0
c
c     pointers for integer work space : iwork
c     ---------------------------------------
      piwo = 1
c     piwo  : pointer for OPMV ( 3*n + neq + 2*nin + 3 )
c
c     Total size of iwork array = ( 3*n + neq + 2*nin + 3 )
c
c     pointers for double precision work space  : dwork
c     -------------------------------------------------
      pdwopt = 1
c     pdwopt : pointer for wopt save vector(n)
      pdwo   = pdwopt + ( n )
c     pdwo   : pointer for internal workspace of OPMV
c              ( n*( n + neq + nin + 9 ) + 2*neq + 4*nin + 4  )
c
c     Total so size of dwork array = n
c                                  + n*(n+neq+nin+9)
c                                  + 2*neq + 4*nin + 4
c         = n*(n + neq + nin + 10) + 2*neq + 4*nin + 4
c
c-----------------------------------------------------------------------
c
c     rdt min/max (rhomin and rhomax)
      CALL EVBOR ( n, rho, rhomin, rhomax )
c
c     research optimal performance subject to Value-at-Risk constraint
c     (upper bound VaR) by dichotomy
c
c     initializations
      vrisav = 0.
      rhopt  = rhomax
      epsvri = 2.*epsdic
      vriopt = -1.0
      infmax = -107
      infqua = 0
      niter  = 0
      xvri   = -1.0
c
c     or number
      ornumber = 2.0/(1.0 + SQRT(5.0))
c
c     loop dichotomy search
      DO WHILE ( ( (niter .LE. maxdic).AND.(epsvri .GT. epsdic) ) .AND.
     &           ( (vrimax .GT. xvri) .OR. (niter .NE. 1) ) )
         niter = niter + 1
c
c        Mean-Variance allocation (Markowitz)
         CALL opmv ( n, cov, rho, rhopt,
     &               neq, nin, ccst, bcst, cinf, csup,
     &               iwork(piwo), dwork(pdwo), wopt, info )
c
c        if wopt(i)< eps then wopt(i) = 0
c        the residual increase max[wopt(i)]
         CALL WZERO ( n, wopt, eps )
c
c        error management
         IF (info .NE. 0) THEN
            IF (info .NE. 1001) THEN
                RETURN
            ELSE
                infqua = 1001
            ENDIF
         ELSE
            infqua = 0
            CALL YV ( n, wopt, dwork(pdwopt) )
         ENDIF
c
c        optimal portfolio ex-ante variance = w'*cov*w
         CALL OVTMCV ( n, cov, wopt, optvar )
c
c        optimal portfolio ex-ante volatility = sqrt(w'*cov*w)
         volat = sqrt(optvar)
c
c        optimal ex-ante mean return = w'*rmean
         CALL XV ( n, wopt, rho, optmean )
c
c        inverse of the N(0,1) cumulative distribution
c        call dinvnr(p,q)
         p  = conflp
         q  = 1. - p
         zc = dinvnr(p,q)
c
c     write p, q, zc
c      open(unit=1,file='zc.txt',status='unknown')
c	write(1,*) p
c	write(1,*) q
c	write(1,*) zc
c	close(unit=1)
c
c        ex-ante Value-at-Risk (upper bound)
c         xvri = - volat * sqrt( conflp/(1.-conflp) )
c         xvri = - volat * sqrt( conflp/(1.-conflp) ) + emret
c
c        ex-ante Normal Value-at-Risk
         xvri = - volat * zc
c         xvri = - volat * zc + emret
c
c     write VaR
c      open(unit=1,file='VaR.txt',status='unknown')
c	write(1,*)  xvri
c	close(unit=1)
c
c        impose that VaR in [-100%, 100%]
         xvri = MIN(xvri, 1.)
         xvri = MAX(xvri, -1.)
c
c        save the last best result
         IF (infqua .EQ. 0) THEN
            vrisav = xvri
         ELSE
            IF (niter .EQ. 1) THEN
               vrisav = xvri
               CALL YV ( n, wopt, dwork(pdwopt) )
               infqua = 0
               info   = 0
            ENDIF
         ENDIF
c
c        stop test
         epsvri = abs(xvri - vriopt)
         vriopt = xvri
c
c        dichotomy (objective return)
         IF (vrimax.lt.vriopt) THEN
            rhomin = rhopt
            infmax = 0
         ELSE
            rhomax = rhopt
         ENDIF
         rhopt = ( rhomax + rhomin ) / 2.
c          rhopt = ornumber * ( rhomax + rhomin )
      ENDDO
c
c     output management
      IF (niter.gt.maxdic) info = 106
      IF (infmax.ne.0) info = -infmax
      IF (infqua.ne.0) THEN
         info = infqua
         CALL YV ( n, dwork(pdwopt), wopt )
         vriopt = vrisav
      ENDIF
c
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
#include "../defines.inc"
!---------------------------------------------------------------*
!
!     Calculate the change in the polymer elastic energy
!     due to the displacement from a MC move
!
subroutine MC_eelas_spider(wlc_p,DEELAS,spider_id,EB,EPAR,EPERP,GAM,ETA)
! values from wlcsim_data
use params, only: wlc_spiders, wlc_U, wlc_R, wlc_UP, wlc_RP

use params, only: dp, wlcsim_params
use MC_wlc, only: E_wlc, E_SSWLC, E_GAUSS
implicit none
type(wlcsim_params), intent(in) :: wlc_p
integer, intent(in) :: spider_id
real(dp), intent(out) :: DEELAS(4)   ! Change in ECOM

!     Polymer properties

real(dp), intent(in) :: EB
real(dp), intent(in) :: EPAR
real(dp), intent(in) :: EPERP
real(dp), intent(in) :: GAM
real(dp), intent(in) :: ETA

integer hip, knee, toe, leg_n, ii, I
integer left_ends(3)  ! left end bead of the three kninks
real(dp) energy_change(4)

! Setup parameters

DEELAS(1) = 0.0_dp ! bending energy
DEELAS(2) = 0.0_dp ! parallel stretch energy
DEELAS(3) = 0.0_dp ! perpendicular stretch energy
DEELAS(4) = 0.0_dp ! WLC_P__TWIST energy

if (WLC_P__RING) then
    print*, "Ring not set up for spider move"
    stop 1
endif

!     Calculate the change in the energy
do leg_n = 1,wlc_spiders(spider_id)%nLegs
    hip = wlc_spiders(spider_id)%legs(1,leg_n)
    knee= wlc_spiders(spider_id)%legs(2,leg_n)
    toe = wlc_spiders(spider_id)%legs(3,leg_n)

    ! note that hip and toe U vectors don't rotate
    ! note also that the knee U vector rotates with the hip
    if (toe>hip) then
        left_ends(1) = hip
        left_ends(2) = knee
        left_ends(3) = toe-1
    else
        left_ends(1) = toe
        left_ends(2) = knee-1
        left_ends(3) = hip-1
    endif
    do ii = 1,3
        I = left_ends(ii)
        ! MC move only affects energy if it's interior to the polymer, since there are only crankshaft moves, and end
        ! crankshafts don't affect polymer
        if (wlc_p%SIMTYPE == 1) then
            print*, "You will need to update this section before use."
            print*, "Finish implementing IT1 and IT2"
            stop 1
            !DEELAS(1) = E_wlc(RP(:,IT1M1), RP(:,IT1), RP(:,IT1P1), EB )
            !DEELAS(1) = DEELAS(1) - E_wlc(R(:,IT1M1), R(:,IT1), R(:,IT1P1), EB)

        elseif (wlc_p%SIMTYPE == 2) then

            !function E_SSWLC(R,RM1,U,UM1,EB,EPAR,EPERP,ETA,GAM)
            energy_change = E_SSWLC(wlc_RP(:,I+1),wlc_RP(:,I),&
                                      wlc_UP(:,I+1),wlc_UP(:,I),&
                                      EB,EPAR,EPERP,ETA,GAM)
            DEELAS = DEELAS + energy_change
            energy_change = E_SSWLC( wlc_R(:,I+1), wlc_R(:,I),&
                                       wlc_U(:,I+1), wlc_U(:,I),&
                                       EB,EPAR,EPERP,ETA,GAM)
            DEELAS = DEELAS - energy_change

        elseif (wlc_p%SIMTYPE == 3) then
            DEELAS(2) = DEELAS(2) + E_GAUSS(wlc_RP(:,I+1),wlc_RP(:,I),EPAR)
            DEELAS(2) = DEELAS(2) - E_GAUSS( wlc_R(:,I+1), wlc_R(:,I),EPAR)
        endif
    enddo
enddo
RETURN
END

!---------------------------------------------------------------*

max_stars_repo_name: SpakowitzLab/BasicWLC

max_issues_repo_name: riscalab/wlcsim

max_forks_repo_name: riscalab/wlcsim


text:
! limiter.f90
! author: sunder

!-----------------------------------------------------------------------
! Routines related to the subcell limiter
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Compute the subcell data lim from a given DG polynomial luh
!-----------------------------------------------------------------------

subroutine get_subcell_data(lim,luh)
    use ader_dg
    implicit none
    ! Argument list
    real :: lim(nVar,nSubLimV(1),nSubLimV(2)), luh(nVar,nDOF(1),nDOF(2))
    ! Local variables
    integer :: iii, jjj
    real :: limy(nVar,nDOF(1),nSubLimV(2))

     ! mapping of uh to the sub-cell limiter along y direction

     limy = 0.0

     do iii = 1, nDOF(1) ! x
         limy(:,iii,:) = matmul( luh(:,iii,:), transpose(uh2lim) )
     end do

     ! mapping of uh to the sub-cell limiter along x direction

     lim = 0.0

     do jjj = 1, nSubLimV(2) ! y
         lim(:,:,jjj) = matmul( limy(:,:,jjj), transpose(uh2lim) )
     end do

 end subroutine get_subcell_data

!-----------------------------------------------------------------------
! Compute the point values lob in the Gauss-Lobatto points from a given DG polynomial luh
!-----------------------------------------------------------------------

subroutine get_lobatto_data(lob,luh)
    use ader_dg
    implicit none
    ! Argument list
    real :: lob(nVar,nDOF(1),nDOF(2)), luh(nVar,nDOF(1),nDOF(2))
    ! Local variables
    integer :: iii, jjj
    real :: loby(nVar,nDOF(1),nDOF(2))

    ! mapping of uh to the Gauss-Lobatto points along y direction

    loby = 0.0

    do iii = 1, nDOF(1) ! x
        loby(:,iii,:) = matmul( luh(:,iii,:), transpose(uh2lob) )
    end do

    ! mapping of uh to the Gauss-Lobatto points along x direction

    lob = 0.0

    do jjj = 1, nDOF(2) ! y
        lob(:,:,jjj) = matmul( loby(:,:,jjj), transpose(uh2lob) )
    end do

end subroutine get_lobatto_data

!-----------------------------------------------------------------------
! Put subcell data back into the solution
!-----------------------------------------------------------------------

subroutine put_subcell_data(luh,lim)
    use ader_dg
    implicit none
    ! Argument list
    real :: lim(nVar,nSubLimV(1),nSubLimV(2)), luh(nVar,nDOF(1),nDOF(2))
    ! Local variables
    integer :: iii, jjj
    real :: luhy(nVar,nSubLimV(1),N+1)

    ! Mapping of the sub-cell limiter to uh along y direction
    luhy = 0.0

    do iii = 1, nSubLimV(1) ! x
        luhy(:,iii,:) = matmul( lim(:,iii,:), transpose(lim2uh) )
    end do

    ! Mapping of the sub-cell limiter to uh along x direction

   luh = 0.0

    do jjj = 1, nDOF(2)  ! y
        luh(:,:,jjj) = matmul( luhy(:,:,jjj), transpose(lim2uh) )
    end do

end subroutine put_subcell_data

!-----------------------------------------------------------------------
! Get the minimum and maximum in the neighbourhood of a cell
!-----------------------------------------------------------------------

subroutine get_min_max
    use ader_dg
    implicit none
    ! Local variables
    integer :: i, j, ii, jj, in, jn,  iVar
    real :: lim(nVar,nSubLimV(1),nSubLimV(2))
    real :: lob(nVar,nDOF(1),nDOF(2))
    real :: Qmin(nVar), Qmax(nVar)

    ! Step 1) Loop over all elements and get the min/max for each variable

    do j = 1, JMAX
        do i = 1, IMAX

            if (Limiter(i,j)%status .eq. 0) then

                call get_subcell_data( lim, uh(:,:,:,i,j) )
                call get_lobatto_data( lob, uh(:,:,:,i,j) )

                do iVar = 1, nVar

                    Limiter(i,j)%lmin_old(iVar) = min( minval(uh(iVar,1:nDOF(1),1:nDOF(2),i,j)), &
                                                   minval(lob(iVar,1:nDOF(1),1:nDOF(2))), &
                                                   minval(lim(iVar,1:nSubLimV(1),1:nSubLimV(2))) )

                    Limiter(i,j)%lmax_old(iVar) = max( maxval(uh(iVar,1:nDOF(1),1:nDOF(2),i,j)), &
                                                   maxval(lob(iVar,1:nDOF(1),1:nDOF(2))), &
                                                   maxval(lim(iVar,1:nSubLimV(1),1:nSubLimV(2))) )

                end do
            else
                lim = Limiter(i,j)%Lh
                do iVar = 1, nVar
                    Limiter(i,j)%lmin_old(iVar) = minval( lim(iVar,1:nSubLimV(1),1:nSubLimV(2)) )
                    Limiter(i,j)%lmax_old(iVar) = maxval( lim(iVar,1:nSubLimV(1),1:nSubLimV(2)) )
                end do
            end if

        end do
    end do

    ! Step 2) Loop over all Voronoi neighbors and get the final min/max for each variable (requires step 1 to be complete)

    do j = 1, JMAX
        do i = 1, IMAX
            Qmin = Limiter(i,j)%lmin_old(:)
            Qmax = Limiter(i,j)%lmax_old(:)

            do jj = -1, 1
                do ii = -1, 1
                    in = neighbor(ii,jj,1,i,j)
                    jn = neighbor(ii,jj,2,i,j)
                    do iVar = 1, nVar
                        Qmin(iVar) = min( Qmin(iVar), Limiter(in,jn)%lmin_old(iVar) )
                        Qmax(iVar) = max( Qmax(iVar), Limiter(in,jn)%lmax_old(iVar) )
                    end do
                end do
            end do

            Limiter(i,j)%lmin(:) = Qmin
            Limiter(i,j)%lmax(:) = Qmax

        end do
    end do

    continue

end subroutine get_min_max

!-----------------------------------------------------------------------
! Check discrete maximum principle (DMP) and physical admissibility
! in the given cell
!-----------------------------------------------------------------------

subroutine DMP(dmpresult,arguh,argLimiter,argmax)
    use ader_dg
    implicit none
    ! Argument list
    logical        :: dmpresult
    real           :: arguh(nVar,nDOF(1),nDOF(2)),argmax
    type(tLimiter) :: argLimiter
    ! Local variables
    integer :: iii, jjj, iVar, iErr
    real    :: lim(nVar,nSubLimV(1),nSubLimV(2))
    real    :: lob(nVar,nDOF(1),nDOF(2))
    real    :: lmin(nVar), lmax(nVar), ldiff, Qout(nVar)

    dmpresult = .true.

    call get_subcell_data( lim, arguh )
    call get_lobatto_data( lob, arguh )

    do iVar = 1, nVar

        if (iVar .gt. 0) then

            lmin(iVar) = min( minval(arguh(iVar,1:nDOF(1),1:nDOF(2))), &
                            minval(lim(iVar,1:nSubLimV(1),1:nSubLimV(2))), &
                            minval(lob(iVar,1:nDOF(1),1:nDOF(2))) )

            lmax(iVar) = max( maxval(arguh(iVar,1:nDOF(1),1:nDOF(2))), &
                            maxval(lim(iVar,1:nSubLimV(1),1:nSubLimV(2))), &
                            maxval(lob(iVar,1:nDOF(1),1:nDOF(2))) )

            ldiff = max( 1.0e-4, 1.0e-3*(argLimiter%lmax(iVar) - argLimiter%lmin(iVar)), argmax )

            if(lmin(iVar) .lt. argLimiter%lmin(iVar)-ldiff) then
                dmpresult = .false.
                return
            end if

            if(lmax(iVar) .gt. argLimiter%lmax(iVar)+ldiff) then
                dmpresult = .false.
                return
            end if

        end if
    end do

    ! Physical detection criteria

    do jjj = 1, nSubLimV(2)
        do iii = 1, nSubLimV(1)
            call PDEAssurePositivity(iErr, lim(:,iii,jjj), Qout )
            if(iErr .lt. 0) then
                dmpresult = .false.
                return
            end if
            !DO iVar = 1, nVar TODO: Add this condition for NaN later
                !IF(ISNAN(lim(iVar,iii,jjj))) THEN
                    !dmpresult = .FALSE.
                    !RETURN
                !ENDIF
            !ENDDO
        end do
    end do

   continue

end subroutine DMP

!-----------------------------------------------------------------------
! Save uh before the limiter is applied
!-----------------------------------------------------------------------

subroutine save_olduh
    use ader_dg
    implicit none
    ! Local variables
    integer  :: i,j

    do j = 1, JMAX
        do i = 1, IMAX
            olduh(:,:,:,i,j) = uh(:,:,:,i,j)
        end do
    end do

    continue

end subroutine save_olduh

!-----------------------------------------------------------------------
! Update the limiter with new status and if required deallocate old
! memory
!-----------------------------------------------------------------------

subroutine update_limiter
    use ader_dg
    implicit none
    ! Local variables
    integer :: i,j

    do j = 1, JMAX
        do i = 1, IMAX

            if(Limiter(i,j)%status .eq. 0) then
                if(Limiter(i,j)%oldstatus .ne. 0) then
                    deallocate( Limiter(i,j)%Lh    )
                    deallocate( Limiter(i,j)%NewLh )
                end if
            else
                Limiter(i,j)%Lh = Limiter(i,j)%NewLh
            end if

            Limiter(i,j)%oldstatus = Limiter(i,j)%status

            continue
        end do
    end do

end subroutine update_limiter

!-----------------------------------------------------------------------
! Allocate memory for the limiter
!-----------------------------------------------------------------------

subroutine allocate_limiter
    use ader_dg
    implicit none
    ! Local variables
    integer    :: i,j

    do j = 1, JMAX
        do i = 1, IMAX
            if( recompute(i,j) .eq. 1) then
                if(Limiter(i,j)%oldstatus .eq. 0) then
                    allocate( Limiter(i,j)%Lh(nVar,nSubLimV(1),nSubLimV(2))    )
                    allocate( Limiter(i,j)%NewLh(nVar,nSubLimV(1),nSubLimV(2)) )
                end if
                Limiter(i,j)%status = 1
            else
                Limiter(i,j)%status = 0
            end if
        end do
    end do

end subroutine allocate_limiter

!-----------------------------------------------------------------------
! Once a cell is marked for recomputation, also mark the face
! neighbours for recomputation
!-----------------------------------------------------------------------

subroutine spread_recompute
    use ader_dg
    implicit none
    ! Local variables
    integer  :: i, j, ii, jj, in, jn
    integer  :: iFace

    do j = 1, JMAX
        do i = 1, IMAX
            if (recompute(i,j) .eq. 1) then
                do iFace = 1, 2*nDim
                    ii = Face2Neigh(1,iFace)
                    jj = Face2Neigh(2,iFace)
                    in = neighbor(ii,jj,1,i,j)
                    jn = neighbor(ii,jj,2,i,j)
                    if (recompute(in,jn) .eq. 0) then
                        recompute(in,jn) = 2
                    end if
                end do
            end if
        end do
    end do

    continue

 end subroutine spread_recompute

 !-----------------------------------------------------------------------
 ! Rusanov flux
 !-----------------------------------------------------------------------

 subroutine RusanovFluxS(F, D, QL, QR, nv)
    USE ader_dg
    implicit none
    ! Argument list
    real, intent(in)  :: QL(nVar), QR(nVar), nv(nDim)
    real, intent(out) :: F(nVar), D(nVar)
    ! Local variables
    real :: smax, LL(nVar), LR(nVar)
    real :: FL(nVar,nDim), FR(nVar,nDim), B(nVar,nVar)

    call PDEEigenvalues(QL, nv, LL)
    call PDEEigenvalues(QR, nv, LR)

    smax = max( maxval(abs(LL)), maxval(abs(LR)) )

    call PDEFlux(QL,FL)
    call PDEFlux(QR,FR)

    F = 0.5*matmul(FL+FR,nv) - 0.5*smax*(QR-QL)

    call RoeMatrix(B, QL, QR, nv)

    D = matmul(B, (QR-QL))

end subroutine RusanovFluxS

!-----------------------------------------------------------------------
! Min-mod limiter
!-----------------------------------------------------------------------

subroutine minmod(c,a,b,nnn)
    implicit none
    ! Argument list
    integer, INTENT(IN)   :: nnn
    real, INTENT(IN)      :: a(nnn), b(nnn)
    real, INTENT(OUT)     :: c(nnn)
    ! Local variables
    integer :: i

    do i = 1, nnn
        if (a(i)*b(i) .le. 0.0) then
            c(i) = 0.0
        else
            if( abs(a(i)) .lt. abs(b(i)) ) then
                c(i) = a(i)
            else
                c(i) = b(i)
            end if
        end if
    end do

end subroutine minmod

!-----------------------------------------------------------------------
! Recompute solution in the subcell
!-----------------------------------------------------------------------

subroutine subcell_recompute
    use ader_dg
    implicit none
    ! Local variables
    integer :: i, j, in, jn, ii, jj
    real    :: ldx(nDim), nv(nDim), nvi(nDim,nDim)
    real    :: subuh0(nVar,(1-nSubLimV(1)):2*nSubLimV(1),(1-nSubLimV(2)):2*nSubLimV(2))
    real    :: subuh1(nVar,nSubLimV(1),nSubLimV(2))
    real    :: slopex(nVar,1-dn(1):nSubLimV(1)+dn(1),1-dn(2):nSubLimV(2)+dn(2))
    real    :: slopey(nVar,1-dn(1):nSubLimV(1)+dn(1),1-dn(2):nSubLimV(2)+dn(2))
    real    :: slopet(nVar,1-dn(1):nSubLimV(1)+dn(1),1-dn(2):nSubLimV(2)+dn(2))
    real    :: wLx(nVar,1-dn(1):nSubLimV(1)+dn(1),1-dn(2):nSubLimV(2)+dn(2))
    real    :: wLy(nVar,1-dn(1):nSubLimV(1)+dn(1),1-dn(2):nSubLimV(2)+dn(2))
    real    :: wRx(nVar,1-dn(1):nSubLimV(1)+dn(1),1-dn(2):nSubLimV(2)+dn(2))
    real    :: wRy(nVar,1-dn(1):nSubLimV(1)+dn(1),1-dn(2):nSubLimV(2)+dn(2))
    real    :: lim(nVar,nSubLimV(1),nSubLimV(2))
    real    :: Qt(nVar), gradQ(nVar,nDim), Src_BgradQ(nVar)
    real    :: FLx(nVar,nDim), FRx(nVar,nDim), FLy(nVar,nDim), FRy(nVar,nDim)
    real    :: Fx(nVar,1:nSubLimV(1)+1,1:nSubLimV(2))
    real    :: Fy(nVar,1:nSubLimV(1),1:nSubLimV(2)+1)
    real    :: fpl(nVar), fmi(nVar), gpl(nVar), gmi(nVar)
    real    :: Dpx(nVar), Dmx(nVar), Dpy(nVar), Dmy(nVar)
    logical :: dmpresult

    nvi = 0.0

    do ii = 1, nDim
        nvi(ii,ii) = 1.0
    end do

    ldx = dx/real(nSubLim)
    Fx = 0.0
    Fy = 0.0

    slopex = 0.0
    slopey = 0.0

    wLx = 0.0
    wRx = 0.0
    wLy = 0.0
    wRy = 0.0

    do j = 1, JMAX
        do i = 1, IMAX
            if (recompute(i,j) > 0) then

                ! Get the initial data, either from a previous limiter stage, or from healthy uh initial data
                ! note that here we only need the values of the face neighbors. For simplicity, we run over all the Voronoi neighbors

                do jj = -dn(2), dn(2)
                    do ii = -dn(1), dn(1)

                        in = neighbor(ii,jj,1,i,j)
                        jn = neighbor(ii,jj,2,i,j)

                        if (Limiter(in,jn)%oldstatus .eq. 0) then
                            call get_subcell_data(lim, olduh(:,:,:,in,jn))
                        else
                            lim = Limiter(in,jn)%Lh
                        end if

                        subuh0(:,1+ii*nSubLimV(1):(ii+1)*nSublimV(1),1+jj*nSubLimV(2):(jj+1)*nSublimV(2)) = lim(:,:,:)

                    end do
                end do

                ! --------------------------------------------------------------
                ! Second order TVD MUSCL reconstruction in space and time
                ! --------------------------------------------------------------

                do jj = 1-dn(2), nSubLimV(2)+dn(2)
                    do ii = 1-dn(1), nSubLimV(1)+dn(1)

                        ! x slopes
                        call minmod( slopex(:,ii,jj), subuh0(:,ii+1,jj)-subuh0(:,ii,jj), subuh0(:,ii,jj)-subuh0(:,ii-1,jj), nVar)
                        ! y slopes
                        call minmod( slopey(:,ii,jj), subuh0(:,ii,jj+1)-subuh0(:,ii,jj), subuh0(:,ii,jj)-subuh0(:,ii,jj-1), nVar)

                        gradQ(:,1) = slopex(:,ii,jj)/ldx(1)
                        gradQ(:,2) = slopey(:,ii,jj)/ldx(2)

                        call PDEFusedSrcNCP(Src_BgradQ, subuh0(:,ii,jj), gradQ)

                        wLx(:,ii,jj) = subuh0(:,ii,jj) - 0.5*slopex(:,ii,jj)
                        wRx(:,ii,jj) = subuh0(:,ii,jj) + 0.5*slopex(:,ii,jj)

                        call PDEFlux(wLx(:,ii,jj), FLx)
                        call PDEFlux(wRx(:,ii,jj), FRx)

                        Qt = Src_BgradQ - (FRx(:,1)-FLx(:,1))/ldx(1)

                        wLy(:,ii,jj) = subuh0(:,ii,jj) - 0.5*slopey(:,ii,jj)
                        wRy(:,ii,jj) = subuh0(:,ii,jj) + 0.5*slopey(:,ii,jj)

                        call PDEFlux(wLy(:,ii,jj), FLy)
                        call PDEFlux(wRy(:,ii,jj), FRy)

                        Qt = Qt - (FRy(:,2)-FLy(:,2))/ldx(2)

                        slopet(:,ii,jj) = Qt

                        wLx(:,ii,jj) = wLx(:,ii,jj) + 0.5*dt*Qt
                        wRx(:,ii,jj) = wRx(:,ii,jj) + 0.5*dt*Qt

                        wLy(:,ii,jj) = wLy(:,ii,jj) + 0.5*dt*Qt
                        wRy(:,ii,jj) = wRy(:,ii,jj) + 0.5*dt*Qt
                    end do
                end do

                do jj = 1, nSubLimV(2)
                    do ii = 1, nSubLimV(1)
                        gradQ(:,1) = slopex(:,ii,jj)/ldx(1)
                        gradQ(:,2) = slopey(:,ii,jj)/ldx(2)

                        call PDEFusedSrcNCP(Src_BgradQ,subuh0(:,ii,jj)+0.5*dt*slopet(:,ii,jj), gradQ)

                        nv = nvi(:,1)

                        call RusanovFluxS(fpl, Dpx, wRx(:,ii,jj), wLx(:,ii+1,jj), nv)
                        call RusanovFluxS(fmi, Dmx, wRx(:,ii-1,jj), wLx(:,ii,jj), nv)

                        nv = nvi(:,2)

                        call RusanovFluxS(gpl, Dpy, wRy(:,ii,jj), wLy(:,ii,jj+1), nv)
                        call RusanovFluxS(gmi, Dmy, wRy(:,ii,jj-1), wLy(:,ii,jj), nv)

                        subuh1(:,ii,jj) = subuh0(:,ii,jj) - dt/ldx(1)*(fpl-fmi) - dt/ldx(2)*(gpl-gmi) &
                                                        - dt/ldx(1)*0.5*(Dpx+Dmx) - dt/ldx(2)*0.5*(Dpy+Dmy) + dt*Src_BgradQ
                    end do
                end do

                if (recompute(i,j) .eq. 1) then ! If a cell is really troubled, then set the limiter status to 1 and set the new subcell data

                    Limiter(i,j)%status = 1

                    ! Save the evolved data in the limiter
                    do jj = 1, nSubLimV(2)
                        do ii = 1, nSubLimV(1)
                            Limiter(i,j)%NewLh(:,ii,jj) = subuh1(:,ii,jj)
                        end do
                    end do

                    ! Put the evolved cell-averaged data back into the DG polynomial
                    call put_subcell_data(uh(:,:,:,i,j),Limiter(i,j)%NewLh)

            else if (recompute(i,j) .eq. 2) then

                    ! Put the evolved cell-averaged data back into the DG polynomial
                    call put_subcell_data(uh(:,:,:,i,j),subuh1)

                    ! If we recompute neighbors of troubled cells, and the resulting reconstructed DG polynomial is troubled, then activate the limiter also in this case
                    call DMP(dmpresult,uh(:,:,:,i,j),Limiter(i,j),0.0)

                    if (dmpresult) then
                        Limiter(i,j)%status = 0
                    else
                        Limiter(i,j)%status = 1
                        if (Limiter(i,j)%oldstatus .eq. 0) then
                        allocate( Limiter(i,j)%Lh(nVar,nSubLimV(1),nSubLimV(2))    )
                        allocate( Limiter(i,j)%NewLh(nVar,nSubLimV(1),nSubLimV(2)) )
                        end if

                        ! Save the evolved data in the limiter

                        do jj = 1, nSubLimV(2)
                            do ii = 1, nSubLimV(1)
                                Limiter(i,j)%NewLh(:,ii,jj) = subuh1(:,ii,jj)
                            end do
                        end do

                    end if
                end if

            end if
        end do
    end do

   continue

end subroutine subcell_recompute

max_stars_repo_name: dasikasunder/ADER-DG

max_issues_repo_name: dasikasunder/ADER-DG

max_forks_repo_name: dasikasunder/ADER-DG


text:
      SUBROUTINE SCALE
*
*
*       Scaling to new units.
*       ---------------------
*
      INCLUDE 'common6.h'
      REAL*8  RSAVE(3,2*KMAX),VSAVE(3,2*KMAX),BSAVE(2*KMAX)
      real*8 G, M_sun, R_sun, pc, Km, Kmps
      real*8 mscale, lscale, vscale
      SAVE RSAVE,VSAVE,BSAVE
*
*       Define physical constants (consistent with IAU & routine UNITS).
      G = 6.6743D-08
      M_sun = 1.9884D+33
      R_sun = 6.955D+10
      pc = 3.0856776D+18
      Km = 1.0D+05
      Kmps = 1.0D+05
*
*       Read virial ratio, rotation scaling factors, tidal radius & SMAX.
      READ (5,*)  Q, VXROT, VZROT, RTIDE, SMAX
*       Note RTIDE should be non-zero for isolated systems (cf. CALL LAGR).
      RSPH2 = RTIDE
      QVIR = Q
*
      ZMASS = 0.0D0
      DO 10 K = 1,3
          CMR(K) = 0.0D0
          CMRDOT(K) = 0.0D0
   10 CONTINUE
*
*       Form total mass and centre of mass displacements.
      DO 30 I = 1,N
          ZMASS = ZMASS + BODY(I)
          DO 25 K = 1,3
              CMR(K) = CMR(K) + BODY(I)*X(K,I)
              CMRDOT(K) = CMRDOT(K) + BODY(I)*XDOT(K,I)
   25     CONTINUE
   30 CONTINUE
*
*       Adjust coordinates and velocities to c.m. rest frame.
      DO 40 I = 1,N
          DO 35 K = 1,3
              X(K,I) = X(K,I) - CMR(K)/ZMASS
              XDOT(K,I) = XDOT(K,I) - CMRDOT(K)/ZMASS
   35     CONTINUE
   40 CONTINUE
*
*       Check optional generation of BH binary/single and other BHs.
      IF (KZ(11).GT.0) THEN
*       Read binary elements and component masses (SEMI = 0 for single BH).
          READ (5,*)  SEMI, ECC, BODY1, BODY2
*       Convert from M_sun to pre-scaled N-body units (ZMASS = N up to now).
          BODY(1) = (BODY1 + BODY2)/ZMBAR
          KSTAR(1) = 14
*       Place a stationary heavy binary/single BH at the centre.
          DO 41 K = 1,3
              X(K,1) = 0.0
              XDOT(K,1) = 0.0
   41     CONTINUE
          IF (SEMI.GT.0.0) THEN
*       Enforce initial KS for massive binary.
              NBIN0 = 1
              NBH0 = 2
              KSTAR(2) = 14
              ILAST = KZ(11) + 1
          ELSE
              NBH0 = 1
              ILAST = KZ(11)
          END IF
*      Add one or more optional free-floating BH.
          IF (KZ(11).GT.1) THEN
              DO 42 I = NBH0+1,ILAST
                  READ (5,*)  BODY(I), (X(K,I),K=1,3), (XDOT(K,I),K=1,3)
                  BODY(I) = BODY(I)/ZMBAR
                  KSTAR(I) = 14
   42         CONTINUE
              NBH0 = ILAST
          END IF
      END IF
*
*       Include procedure for primordial binaries & singles from fort.10.
      IF ((KZ(22).EQ.4.OR.KZ(22).EQ.-1).AND.(NBIN0.GT.0)) THEN
*       Save the binaries for re-creating two-body motion after scaling.
          DO 45 I = 1,2*NBIN0
              BSAVE(I) = BODY(I)
              DO 44 K = 1,3
                  RSAVE(K,I) = X(K,I)
                  VSAVE(K,I) = XDOT(K,I)
   44         CONTINUE
   45     CONTINUE
*       Form the c.m. of each binary for standard scaling to E = -0.25.
          DO 47 I = 1,NBIN0
              I1 = 2*I - 1
              I2 = 2*I
*       Prevent over-writing BODY(I) when I = 1 and 2*I - 1 = 1.
              B1 = BODY(I1)
              B2 = BODY(I2)
              ZMB  = B1 + B2
              BODY(I) = ZMB
              DO 46 K = 1,3
                  X(K,I) = (B1*X(K,I1) + B2*X(K,I2))/ZMB
                  XDOT(K,I) = (B1*XDOT(K,I1) +
     &                         B2*XDOT(K,I2))/ZMB
   46         CONTINUE
   47     CONTINUE
*
*       Move any single stars up for scaling together with c.m. particles.
          NSING = N - 2*NBIN0
*       Note assumption ZMASS = 1 including any singles (skips DO 50 loop).
          DO 49 I = 1,NSING
             BODY(NBIN0+I) = BODY(2*NBIN0+I)
             DO 48 K = 1,3
                 X(K,NBIN0+I) = X(K,2*NBIN0+I)
                 XDOT(K,NBIN0+I) = XDOT(K,2*NBIN0+I)
   48        CONTINUE
   49     CONTINUE
*       Re-define N & NTOT for total energy calculation (NTOT = N for ZKIN).
          N = NBIN0 + NSING
          NTOT = N
*
          write (6,491)  NBIN0
  491     format (/,12X,'SCALE: ', I5,
     &                  ' Binaries converted to barycentres')
      END IF
*       Save total number of C.Ms.
      NCM = NTOT
*
*       Skip scaling of masses for unscaled upload or planetesimal disk.
      IF (KZ(22).GT.2.OR.KZ(22).EQ.-1.OR.KZ(5).EQ.3) GO TO 52
*
*       Scale masses to standard units of <M> = 1/N and set total mass.
      DO 50 I = 1,N
          BODY(I) = BODY(I)/ZMASS
   50 CONTINUE
      ZMASS = 1.0
*
*       Obtain the total kinetic & potential energy.
   52 CALL ENERGY
*
*       Check option for astrophysical units.
      if (KZ(22).EQ.-1) then
*       Save K.E. and P.E. of C.Ms for TCR & TRH (units of M_sun, pc & km/s).
          ZCM = ZKIN
          PCM = POT
*
*       Evaluate total energy in physical units (M_sun pc^2 s^-2).
*       G (g^-1 cm^3 s^-2) ==> G (M_sun^-1 pc^3 s^-2)
          G = G/(pc**3/M_sun)
*       Form kinetic energy as KE (M_sun Km^2 s^-2) ==> KE (M_sun pc^2 s^-2).
          ZKIN = ZKIN*Km**2/pc**2
*       PE term = (M_sun^-1 pc^3 s^-2)*(M_sun^2 pc^-1) = (M_sun pc^2 s^-2).
          EPH = ZKIN - G*POT
          IF (EPH.GT.0.0) EPH = -EPH
*       Set conversion factors (N-body to M_sun, pc, km/s; Heggie & Mathieu).
          mscale = ZMASS
          lscale = G*mscale**2*(-0.25/EPH)
          vscale = sqrt(G*mscale/lscale)*(pc/Km)
*       Specify virial radius and mean mass (pc & M_sun).
          RBAR = lscale
          ZMBAR = ZMASS/(N+NBIN0)
*
          write (6,53)  mscale, lscale, vscale
   53     format (/,12X,'SCALE: Conversion factors: MSCALE = ',F10.3,
     &                  ' M_sun;  LSCALE = ',F6.3,' pc;  VSCALE = ',
     &                    F6.3,' km/sec')
      end if
*
*       Use generalized virial theorem for external tidal field.
      IF (KZ(14).GT.0) THEN
          AZ = 0.0D0
          DO 55 I = 1,N
              AZ = AZ + BODY(I)*(X(1,I)*XDOT(2,I) - X(2,I)*XDOT(1,I))
   55     CONTINUE
          IF (KZ(14).EQ.1) THEN
*       Use Chandrasekhar eq. (5.535) for virial ratio (rotating frame only).
              VIR = POT - 2.0*(ETIDE + 0.5*TIDAL(4)*AZ)
          ELSE
              VIR = POT - 2.0*ETIDE
          END IF
      ELSE
          VIR = POT
      END IF
*
*       Allow two optional ways of skipping standard velocity scaling.
      IF (KZ(22).EQ.3.OR.KZ(5).EQ.2.OR.KZ(5).EQ.3) THEN
          QV = SQRT(Q*VIR/ZKIN)
          E0 = ZKIN*QV**2 - POT + ETIDE
          SX = 1.0
*       Rescale velocities to new masses for two Plummer spheres.
          IF (KZ(5).EQ.2) THEN
              ZKIN = 0.0
              DO 57 I = 1,N
                  DO 56 K = 1,3
                      XDOT(K,I) = XDOT(K,I)*QV
                      ZKIN = ZKIN + 0.5*BODY(I)*XDOT(K,I)**2
   56             CONTINUE
   57         CONTINUE
              E0 = ZKIN - POT + ETIDE
              Q = ZKIN/POT
              WRITE (6,59)  E0, ZKIN/POT
   59         FORMAT (/,12X,'UNSCALED ENERGY    E =',F10.6,
     &                                       '  Q =',F6.2)
          ELSE
              IF (KZ(5).EQ.3) E0 = ZKIN - POT
              WRITE (6,54)  E0
   54         FORMAT (/,12X,'UNSCALED ENERGY    E =',F10.6)
          END IF
      ELSE IF(KZ(22).NE.-1) THEN
*       Scale non-zero velocities by virial theorem ratio.
          IF (ZKIN.GT.0.0D0) THEN
              QV = SQRT(Q*VIR/ZKIN)
              DO 60 I = 1,N
                  DO 58 K = 1,3
                      XDOT(K,I) = XDOT(K,I)*QV
   58             CONTINUE
   60         CONTINUE
          ELSE
              QV = 1.0
          END IF
*
*       Scale total energy to standard units (E = -0.25 for Q < 1).
          E0 = -0.25
*       Include case of hot system inside reflecting boundary.
          IF (KZ(29).GT.0.AND.Q.GT.1.0) E0 = 0.25
*         ETOT = (Q - 1.0)*POT
          ETOT = ZKIN*QV**2 - POT + ETIDE
*       Note that final ETOT will differ from -0.25 since ETIDE = 0.
          IF (Q.LT.1.0) THEN
              SX = E0/ETOT
          ELSE
              SX = 1.0
          END IF
*
          WRITE (6,65)  SX, ETOT, BODY(1), BODY(N), ZMASS/FLOAT(N)
   65     FORMAT (/,12X,'SCALING:    SX =',F6.2,'  E =',1PE10.2,
     &                  '  M(1) =',E9.2,'  M(N) =',E9.2,'  <M> =',E9.2)
*
*       Scale coordinates & velocities to the new units.
          DO 70 I = 1,N
              DO 68 K = 1,3
                  X(K,I) = X(K,I)/SX
                  XDOT(K,I) = XDOT(K,I)*SQRT(SX)
   68         CONTINUE
   70     CONTINUE
*
*       Set current energies to consistent values (may be needed in XTRNL0).
          ZKIN = ZKIN*QV**2*SX
          POT = POT*SX
      END IF
*
*       Perform second stage of the optional binary uploading procedure.
      IF ((KZ(22).EQ.4.OR.KZ(22).EQ.-1).AND.(NBIN0.GT.0)) THEN
*       Place any singles last and re-create the original binaries.
          DO 72 I = NSING,1,-1
              L1 = NBIN0
              L2 = 2*NBIN0
              BODY(L2+I) = BODY(L1+I)
              DO 71 K = 1,3
                  X(K,L2+I) = X(K,L1+I)
                  XDOT(K,L2+I) = XDOT(K,L1+I)
   71         CONTINUE
   72     CONTINUE
*       Split each c.m. into binary components using saved quantities.
          DO 74 I = 1,NBIN0
              IP = NBIN0 - I + 1
              BODY(2*IP-1) = BSAVE(2*IP-1)
              BODY(2*IP) = BSAVE(2*IP)
              ZMB = BODY(2*IP-1) + BODY(2*IP)
*       Use the original relative motion for unscaled two-body elements.
              DO 73 K = 1,3
                  XREL = RSAVE(K,2*IP-1) - RSAVE(K,2*IP)
                  VREL = VSAVE(K,2*IP-1) - VSAVE(K,2*IP)
*       Prevent over-writing of first location when IP = 1 and 2*IP-1 = 1.
                  X1 = X(K,IP)
                  V1 = XDOT(K,IP)
                  X(K,2*IP-1) = X(K,IP) + BSAVE(2*IP)*XREL/ZMB
                  X(K,2*IP) = X1 - BSAVE(2*IP-1)*XREL/ZMB
                  XDOT(K,2*IP-1) = XDOT(K,IP) + BSAVE(2*IP)*VREL/ZMB
                  XDOT(K,2*IP) = V1 - BSAVE(2*IP-1)*VREL/ZMB
   73         CONTINUE
   74     CONTINUE
*       Set final values of the particle number (increase from N + NBIN0).
          N = N + NBIN0
          NZERO = N
          NTOT = N
*
          write (6,741)  NBIN0
  741     format (/,12X,'SCALE: ', I5, ' Binaries restored')
      END IF
*
*       Scale masses, coordinates and velocities for input in physical units.
      if (KZ(22).EQ.-1) then
*       Re-calculate the scaled total mass.
          ZMASS = 0.0D0
          DO 743 I = 1,NTOT
             BODY(I) = BODY(I)/mscale
             ZMASS = ZMASS + BODY(I)
             DO 742 K = 1,3
                X(K,I) = X(K,I)/lscale
                XDOT(K,I) = XDOT(K,I)/vscale
  742        continue
  743     continue
*
*       Perform energy check after scaling.
          CALL ENERGY
          write (6,745)  ZKIN - POT
  745     format (/,12X,'SCALE: Scaled total energy: ',F8.5)
*       Set energy of singles & C.Ms only in N-body units for TCR & TRH.
*       escale = mscale*(vscale**2)
          SX = 1.0
          ZKIN = ZCM/(mscale*vscale**2)
          POT = PCM*lscale/mscale**2
          E0 = ZKIN - POT
          write (6,746)  E0
  746     format (/,12X,'SCALE: Scaled C.M. energy: ',F8.5)
      end if
*
*       Introduce optional BH binary (SEMI in scaled N-body units).
      IF (KZ(11).GT.0) THEN
          IF (SEMI.GT.0.0) THEN
              ZMB = BODY(1)
*       Split c.m. into components with given mass ratio.
              BODY(1) = BODY1/(BODY1 + BODY2)*ZMB
              BODY(2) = BODY2/(BODY1 + BODY2)*ZMB
              RAP = SEMI*(1.0 + ECC)
              VAP = SQRT(ZMB*(1.0 - ECC)/(SEMI*(1.0 + ECC)))
*       Specify two-body motion for SEMI & ECC in x-y plane.
              X(1,2) = X(1,1) - BODY(1)*RAP/ZMB
              X(1,1) = X(1,1) + BODY(2)*RAP/ZMB
              XDOT(2,2) = XDOT(2,1) - BODY(1)*VAP/ZMB
              XDOT(2,1) = XDOT(2,1) + BODY(2)*VAP/ZMB
              DO 75 I = 1,2
                  X(2,I) = 0.0
                  X(3,I) = 0.0
                  XDOT(1,I) = 0.0
                  XDOT(3,I) = 0.0
   75         CONTINUE
          END IF
          DO 77 I = 1,ILAST
              WRITE (6,76)  BODY(I), (X(K,I),K=1,3), (XDOT(K,I),K=1,3)
   76         FORMAT (12X,'BH COMPONENTS    M X XDOT ',
     &                                      1P,E10.2,0P,2(2X,3F7.3))
   77     CONTINUE
      END IF
*
*       Check whether to include rotation (VXROT = 0 in standard case).
      IF (VXROT.GT.0.0D0) THEN
*
*       Set angular velocity for retrograde motion (i.e. star clusters).
          OMEGA = -SX*SQRT(ZMASS*SX)
          WRITE (6,78)  VXROT, VZROT, OMEGA
   78     FORMAT (/,12X,'VXROT =',F6.2,'  VZROT =',F6.2,
     &                                                 '  OMEGA =',F7.2)
*
*       Add solid-body rotation about Z-axis (reduce random velocities).
          DO 80 I = 1,N
              XDOT(1,I) = XDOT(1,I)*VXROT - X(2,I)*OMEGA
              XDOT(2,I) = XDOT(2,I)*VXROT + X(1,I)*OMEGA
              XDOT(3,I) = XDOT(3,I)*VZROT
   80     CONTINUE
      END IF
*
*       Check option for writing the initial conditions on unit #10.
      IF (KZ(22).EQ.1) THEN
          DO 85 I = 1,N
              WRITE (10,84)  BODY(I), (X(K,I),K=1,3), (XDOT(K,I),K=1,3)
   84         FORMAT (1P,7E14.6)
   85     CONTINUE
      END IF
*
*       Check for writing initial conditions (physical units) on unit #99.
      IF (KZ(22).EQ.-1.OR.KZ(22).EQ.5) THEN
          DO 86 I = 1,N
              write (99,84)  BODY(I)*mscale, (X(K,I)*lscale,K=1,3),
     &                       (XDOT(K,I)*vscale,K=1,3)
              BODY0(I) = BODY(I)
   86     CONTINUE
      END IF
*
*       Check option for reading initial subsystems (solar masses).
      IF (KZ(24).GT.0) THEN
          K = KZ(24)
          DO 90 I = 1,K
              READ (5,*)  BODY(I), (X(J,I),J=1,3), (XDOT(J,I),J=1,3)
              BODY(I) = BODY(I)/(ZMBAR*FLOAT(N))
   90     CONTINUE
      END IF
*
*       Set initial crossing time in scaled units.
      TCR = ZMASS**2.5/(2.0D0*ABS(E0))**1.5
      TCR0 = TCR
*
*       Obtain approximate half-mass radius after scaling.
      RSCALE = 0.5*ZMASS**2/POT
*       Set square radius of reflecting sphere (used with option 29).
      RSPH2 = (RSPH2*RSCALE)**2
*       Form equilibrium rms velocity (temporarily defined as VC).
      VC = SQRT(2.0D0*ABS(E0)/ZMASS)
*
*       Check for general binary search of initial condition.
*     IF (KZ(4).GT.0) THEN
*         CALL EVOLVE(0,0)
*     END IF
*
*       Print half-mass relaxation time & equilibrium crossing time.
      A1 = FLOAT(NCM)
      TRH = 4.0*TWOPI/3.0*(VC*RSCALE)**3/(15.4*ZMASS**2*LOG(A1)/A1)
      WRITE (6,95)  TRH, TCR, 2.0*RSCALE/VC, SMAX
   95 FORMAT (/,12X,'TIME SCALES:   TRH =',1P,E8.1,'  TCR =',0P,F6.2,
     &                           '  2<R>/<V> =',F6.2,'  SMAX =',F7.3,/)
*
      RETURN
*
      END


text:
*DECK XERHLT
      SUBROUTINE XERHLT (MESSG)
C***BEGIN PROLOGUE  XERHLT
C***SUBSIDIARY
C***PURPOSE  Abort program execution and print error message.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3C
C***TYPE      ALL (XERHLT-A)
C***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        ***Note*** machine dependent routine
C        XERHLT aborts the execution of the program.
C        The error message causing the abort is given in the calling
C        sequence, in case one needs it for printing on a dayfile,
C        for example.
C
C     Description of Parameters
C        MESSG is as in XERMSG.
C
C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
C                 Error-handling Package, SAND82-0800, Sandia
C                 Laboratories, 1982.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900206  Routine changed from user-callable to subsidiary.  (WRB)
C   900510  Changed calling sequence to delete length of character
C           and changed routine name from XERABT to XERHLT.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XERHLT
      CHARACTER*(*) MESSG
C***FIRST EXECUTABLE STATEMENT  XERHLT
      write(*,*) 'stopping...'
      STOP
      END


! { dg-do run }

  use omp_lib

  double precision :: d, e
  logical :: l
  integer (kind = omp_lock_kind) :: lck
  integer (kind = omp_nest_lock_kind) :: nlck

  d = omp_get_wtime ()

  call omp_init_lock (lck)
  call omp_set_lock (lck)
  if (omp_test_lock (lck)) call abort
  call omp_unset_lock (lck)
  if (.not. omp_test_lock (lck)) call abort
  if (omp_test_lock (lck)) call abort
  call omp_unset_lock (lck)
  call omp_destroy_lock (lck)

  call omp_init_nest_lock (nlck)
  if (omp_test_nest_lock (nlck) .ne. 1) call abort
  call omp_set_nest_lock (nlck)
  if (omp_test_nest_lock (nlck) .ne. 3) call abort
  call omp_unset_nest_lock (nlck)
  call omp_unset_nest_lock (nlck)
  if (omp_test_nest_lock (nlck) .ne. 2) call abort
  call omp_unset_nest_lock (nlck)
  call omp_unset_nest_lock (nlck)
  call omp_destroy_nest_lock (nlck)

  call omp_set_dynamic (.true.)
  if (.not. omp_get_dynamic ()) call abort
  call omp_set_dynamic (.false.)
  if (omp_get_dynamic ()) call abort

  call omp_set_nested (.true.)
  if (.not. omp_get_nested ()) call abort
  call omp_set_nested (.false.)
  if (omp_get_nested ()) call abort

  call omp_set_num_threads (5)
  if (omp_get_num_threads () .ne. 1) call abort
  if (omp_get_max_threads () .ne. 5) call abort
  if (omp_get_thread_num () .ne. 0) call abort
  call omp_set_num_threads (3)
  if (omp_get_num_threads () .ne. 1) call abort
  if (omp_get_max_threads () .ne. 3) call abort
  if (omp_get_thread_num () .ne. 0) call abort
  l = .false.
!$omp parallel reduction (.or.:l)
  l = omp_get_num_threads () .ne. 3
  l = l .or. (omp_get_thread_num () .lt. 0)
  l = l .or. (omp_get_thread_num () .ge. 3)
!$omp master
  l = l .or. (omp_get_thread_num () .ne. 0)
!$omp end master
!$omp end parallel
  if (l) call abort

  if (omp_get_num_procs () .le. 0) call abort
  if (omp_in_parallel ()) call abort
!$omp parallel reduction (.or.:l)
  l = .not. omp_in_parallel ()
!$omp end parallel
!$omp parallel reduction (.or.:l) if (.true.)
  l = .not. omp_in_parallel ()
!$omp end parallel

  e = omp_get_wtime ()
  if (d .gt. e) call abort
  d = omp_get_wtick ()
  ! Negative precision is definitely wrong,
  ! bigger than 1s clock resolution is also strange
  if (d .le. 0 .or. d .gt. 1.) call abort
end


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! This is for getting & storing the alpha and kappa tables

module akmod
  implicit none

  double precision :: muvect(100), alphavect(100), kappavect(100)
  integer, parameter :: nlinesinaktable =100 !because lazy
contains

  subroutine interp1(xin,yin,lenxin,xout,yout)
    !1d interpolation assuming monotonically increasing vector
    integer, intent(in) :: lenxin
    double precision, intent(in) :: xin(lenxin), yin(lenxin), xout
    double precision :: yout
    integer :: i

    if (xout .lt. xin(1)) then
      print*,"xout ", xout, "xin min ", xin(1)
      stop "Error in interpolation: xout < min(xin)"
    end if
    if (xout .gt. xin(lenxin)) then
    print*,"xout ", xout, "xin max ", xin(lenxin), lenxin
    stop "Error in interpolation: xout > max(xin)"
  end if

    i = 1
    do while (xout .gt. xin(i))
      i = i+1
    end do

    yout = yin(i-1)+(yin(i)-yin(i-1))/(xin(i)-xin(i-1))*(xout-xin(i-1))
    return
  end subroutine interp1

end module akmod

subroutine get_alpha_kappa(nq,nv)
  use akmod
  implicit none
  !v^2nv, q^2nq
  integer, intent(in) :: nq, nv
  integer i
  character*300 afilename, kfilename

  open(99,file="ak_files/mVect.dat")
  do i=1,nlinesinaktable
  read(99,*) muVect(i)
  end do
  close(99)

  if ((nq .eq. 0) .and. (nv .eq. 0)) then
     afilename = "ak_files/aVect_0.dat"
     kfilename = "ak_files/kVect_0.dat"
   else if (nq == 1) then
     afilename = "ak_files/aVect_q2.dat"
     kfilename = "ak_files/kVect_q2.dat"
   else if (nq == 2) then
     afilename = "ak_files/aVect_q4.dat"
     kfilename = "ak_files/kVect_q4.dat"
   else if (nq == -1) then
     afilename = "ak_files/aVect_qm2.dat"
     kfilename = "ak_files/kVect_qm2.dat"
   else if (nv == 1) then
     afilename = "ak_files/aVect_v2.dat"
     kfilename = "ak_files/kVect_v2.dat"
   else if (nv == 2) then
     afilename = "ak_files/aVect_v4.dat"
     kfilename = "ak_files/kVect_v4.dat"
   else if (nv == -1) then
     afilename = "ak_files/aVect_vm2.dat"
     kfilename = "ak_files/kVect_vm2.dat"
   end if

  open(99,file=afilename)
  open(95,file=kfilename)
  do i=1,nlinesinaktable
  read(99,*) alphavect(i)
  read(95,*) kappavect(i)
  end do
  close(99)
  close(95)

  return
end subroutine get_alpha_kappa