+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (itype(i-2,1).eq.0) then
+ iti = nloctyp
+ else
+ iti = itype2loc(itype(i-2,1))
+ endif
+ else
+ iti=nloctyp
+ endif
+!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ iti1 = itype2loc(itype(i-1,1))
+ else
+ iti1=nloctyp
+ endif
+! print *,i,itype(i-2,1),iti
+#ifdef NEWCORR
+ cost1=dcos(theta(i-1))
+ sint1=dsin(theta(i-1))
+ sint1sq=sint1*sint1
+ sint1cub=sint1sq*sint1
+ sint1cost1=2*sint1*cost1
+! print *,"cost1",cost1,theta(i-1)
+!c write (iout,*) "bnew1",i,iti
+!c write (iout,*) (bnew1(k,1,iti),k=1,3)
+!c write (iout,*) (bnew1(k,2,iti),k=1,3)
+!c write (iout,*) "bnew2",i,iti
+!c write (iout,*) (bnew2(k,1,iti),k=1,3)
+!c write (iout,*) (bnew2(k,2,iti),k=1,3)
+ k=1
+! print *,bnew1(1,k,iti),"bnew1"
+ do k=1,2
+ b1k=bnew1(1,k,iti)+(bnew1(2,k,iti)+bnew1(3,k,iti)*cost1)*cost1
+! print *,b1k
+! write(*,*) shape(b1)
+! if(.not.allocated(b1)) print *, "WTF?"
+ b1(k,i-2)=sint1*b1k
+!
+! print *,b1(k,i-2)
+
+ gtb1(k,i-2)=cost1*b1k-sint1sq*&
+ (bnew1(2,k,iti)+2*bnew1(3,k,iti)*cost1)
+! print *,gtb1(k,i-2)
+
+ b2k=bnew2(1,k,iti)+(bnew2(2,k,iti)+bnew2(3,k,iti)*cost1)*cost1
+ b2(k,i-2)=sint1*b2k
+! print *,b2(k,i-2)
+
+ gtb2(k,i-2)=cost1*b2k-sint1sq*&
+ (bnew2(2,k,iti)+2*bnew2(3,k,iti)*cost1)
+! print *,gtb2(k,i-2)
+
+ enddo
+! print *,b1k,b2k
+ do k=1,2
+ aux=ccnew(1,k,iti)+(ccnew(2,k,iti)+ccnew(3,k,iti)*cost1)*cost1
+ cc(1,k,i-2)=sint1sq*aux
+ gtcc(1,k,i-2)=sint1cost1*aux-sint1cub*&
+ (ccnew(2,k,iti)+2*ccnew(3,k,iti)*cost1)
+ aux=ddnew(1,k,iti)+(ddnew(2,k,iti)+ddnew(3,k,iti)*cost1)*cost1
+ dd(1,k,i-2)=sint1sq*aux
+ gtdd(1,k,i-2)=sint1cost1*aux-sint1cub*&
+ (ddnew(2,k,iti)+2*ddnew(3,k,iti)*cost1)
+ enddo
+! print *,"after cc"
+ cc(2,1,i-2)=cc(1,2,i-2)
+ cc(2,2,i-2)=-cc(1,1,i-2)
+ gtcc(2,1,i-2)=gtcc(1,2,i-2)
+ gtcc(2,2,i-2)=-gtcc(1,1,i-2)
+ dd(2,1,i-2)=dd(1,2,i-2)
+ dd(2,2,i-2)=-dd(1,1,i-2)
+ gtdd(2,1,i-2)=gtdd(1,2,i-2)
+ gtdd(2,2,i-2)=-gtdd(1,1,i-2)
+! print *,"after dd"
+
+ do k=1,2
+ do l=1,2
+ aux=eenew(1,l,k,iti)+eenew(2,l,k,iti)*cost1
+ EE(l,k,i-2)=sint1sq*aux
+ gtEE(l,k,i-2)=sint1cost1*aux-sint1cub*eenew(2,l,k,iti)
+ enddo
+ enddo
+ EE(1,1,i-2)=EE(1,1,i-2)+e0new(1,iti)*cost1
+ EE(1,2,i-2)=EE(1,2,i-2)+e0new(2,iti)+e0new(3,iti)*cost1
+ EE(2,1,i-2)=EE(2,1,i-2)+e0new(2,iti)*cost1+e0new(3,iti)
+ EE(2,2,i-2)=EE(2,2,i-2)-e0new(1,iti)
+ gtEE(1,1,i-2)=gtEE(1,1,i-2)-e0new(1,iti)*sint1
+ gtEE(1,2,i-2)=gtEE(1,2,i-2)-e0new(3,iti)*sint1
+ gtEE(2,1,i-2)=gtEE(2,1,i-2)-e0new(2,iti)*sint1
+! print *,"after ee"
+
+!c b1tilde(1,i-2)=b1(1,i-2)
+!c b1tilde(2,i-2)=-b1(2,i-2)
+!c b2tilde(1,i-2)=b2(1,i-2)
+!c b2tilde(2,i-2)=-b2(2,i-2)
+#ifdef DEBUG
+ write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
+ write(iout,*) 'b1=',(b1(k,i-2),k=1,2)
+ write(iout,*) 'b2=',(b2(k,i-2),k=1,2)
+ write (iout,*) 'theta=', theta(i-1)
+#endif
+#else
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ iti = itype2loc(itype(i-2,1))
+ else
+ iti=nloctyp
+ endif
+!c write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti
+!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ iti1 = itype2loc(itype(i-1,1))
+ else
+ iti1=nloctyp
+ endif
+! print *,i,iti
+ b1(1,i-2)=b(3,iti)
+ b1(2,i-2)=b(5,iti)
+ b2(1,i-2)=b(2,iti)
+ b2(2,i-2)=b(4,iti)
+ do k=1,2
+ do l=1,2
+ CC(k,l,i-2)=ccold(k,l,iti)
+ DD(k,l,i-2)=ddold(k,l,iti)
+ EE(k,l,i-2)=eeold(k,l,iti)
+ enddo
+ enddo
+#endif
+ b1tilde(1,i-2)= b1(1,i-2)
+ b1tilde(2,i-2)=-b1(2,i-2)
+ b2tilde(1,i-2)= b2(1,i-2)
+ b2tilde(2,i-2)=-b2(2,i-2)
+!c
+ Ctilde(1,1,i-2)= CC(1,1,i-2)
+ Ctilde(1,2,i-2)= CC(1,2,i-2)
+ Ctilde(2,1,i-2)=-CC(2,1,i-2)
+ Ctilde(2,2,i-2)=-CC(2,2,i-2)
+!c
+ Dtilde(1,1,i-2)= DD(1,1,i-2)
+ Dtilde(1,2,i-2)= DD(1,2,i-2)
+ Dtilde(2,1,i-2)=-DD(2,1,i-2)
+ Dtilde(2,2,i-2)=-DD(2,2,i-2)
+ enddo
+#ifdef PARMAT
+ do i=ivec_start+2,ivec_end+2
+#else
+ do i=3,nres+1
+#endif
+