c set of subroutines for the cz*.for programs:
!	subroutine unit(a,n)
!       subroutine vector(v,cte,n)
!	subroutine diagop(a,v,n)
!	subroutine diago(a,v,n)             diagonal matrix with v
!	subroutine dyago(a,v,n)             diagonal matrix with inverse of v
!	subroutine invdiag(a,b,n)           inverse of diagonal matrix 
!	subroutine sypvvv(a,b,c,d,n)        suma y prod de 3 vectores, muy comun
!	subroutine sypvmv(v,w,b,u,n)        suma y prod de 3 vectores, muy comun
!	subroutine mulmvv(w,b,u,v,n)        prod matriz vector vector
!	subroutine muymvv(w,b,u,v,n)        prod matriz (inv.vector) vector
!	subroutine samem (a,m,n)
!	subroutine samemcore (a,m,n,n-2)    extract core of matrix
!	subroutine samemsp (a,m,n)
!	subroutine samevsp (v,w,n)
!	subroutine samev (v,w,n)
!	subroutine samevcore (v,w,n,n-2)    extract core of vector
! no	subroutine operaux (a,n,  b,c,d,e, ll,mm,dd,maux1,maux2) 
! no	subroutine invmcore (a,acore,n,  dd,ll,mm) 
!	subroutine mulmv(a,b,c,n)
!	subroutine mulvmv(a,u,b,v,n)
!	subroutine mulmm(a,b,c,n)
!	subroutine summm(a,b,c,n)
!	subroutine resmm(a,b,c,n)
!	subroutine mulvv(a,b,c,n)
!	subroutine sumvv(a,b,c,n)
!	subroutine sumvvv(a,b,c,d,n)
!	subroutine resvv(a,b,c,n)
!	subroutine zerom(a,n)
!	subroutine zeromsp (a,n)
!	subroutine zero4m(a,b,c,d,n)
!	subroutine zero4msp(a,b,c,d,n)
!	subroutine zero3m(a,b,c,n)
!	subroutine zero3msp(a,b,c,n)
!	subroutine zero2m(a,b,n)
!	subroutine zero2msp(a,b,n)
!	subroutine zerov(a,n)
!	subroutine zerovdim3(a,n1,n2,n3)  ! sustituye a zerojt de cristina
!	subroutine zero4v(a,b,c,d,n)
!	subroutine zero3v(a,b,c,n)
!	subroutine zero2v(a,b,n)
!	subroutine zerovsp(a,n)
!	subroutine zero4vsp(a,b,c,d,n)
!	subroutine zero3vsp(a,b,c,n)
!	subroutine zero2vsp(a,b,n)
!
!
!
!   May-05 Sustituimos todos los zerojt de cristina por las subrutinas
!          genericas zerov***
!
c	***********************************************************************
	subroutine unit(a,n)
c	store the unit value in the diagonal of a 
c	***********************************************************************
	real*8 a(n,n)
	integer n,i,j,k
	do 1,i=2,n-1
	  do 2,j=2,n-1
	    if(i.eq.j) then
	      a(i,j) = 1.d0
	    else
	      a(i,j)=0.0d0
	    end if
 2	  continue
 1	continue
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do
	return
	end

!       subroutine vector(v,cte,n)
c	***********************************************************************
	subroutine vector(v,cte,n)
c	build a vector by storing the value cte in all its elements
c	***********************************************************************
	real*8 v(n),cte
	integer n,i
	do 1,i=1,n
	  v(i) = cte
 1	continue
	return
	end
c	***********************************************************************
	subroutine diagop(a,v,n)
c	store the core of v in the diagonal elements of the square matrix a
c	***********************************************************************
	real*8 a(n,n),v(n+2)
	integer n,i,j,k
	do 1,i=2,n-1
	  do 2,j=2,n-1
	    if(i.eq.j) then
	      a(i,j) = v(i+1)
	    else
	      a(i,j)=0.0d0
	    end if
 2	  continue
 1	continue
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do
	return
	end
c	***********************************************************************
	subroutine diago(a,v,n)
c	store the vector v in the diagonal elements of the square matrix a
c	***********************************************************************
	implicit none

	integer n,i,j,k
	real*8 a(n,n),v(n)

	do 1,i=2,n-1
	  do 2,j=2,n-1
	    if(i.eq.j) then
	      a(i,j) = v(i)
	    else
	      a(i,j)=0.0d0
	    end if
 2	  continue
 1	continue
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do
	return
	end  
c	***********************************************************************
	subroutine dyago(a,v,n)
c	store the inverse of v in the diagonal elements of the square matrix a
c	***********************************************************************
	implicit none

	integer n,i,j,k
	real*8 a(n,n),v(n)

	do 1,i=2,n-1
	  do 2,j=2,n-1
	    if(i.eq.j) then
	      a(i,j) = 1.d0/v(i)
	    else
	      a(i,j)=0.0d0
	    end if
 2	  continue
 1	continue
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do
	return
	end

c	***********************************************************************
	subroutine samem (a,m,n)
c	store the matrix m in the matrix a 
c	***********************************************************************
	real*8 a(n,n),m(n,n)
	integer n,i,j,k
	do 1,i=2,n-1
	  do 2,j=2,n-1
	      a(i,j) = m(i,j)  
2	  continue
1	continue 	
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do
	return
	end 
c	***********************************************************************
	subroutine samemcore (a,b,m,n)
c	store the matrix m in the matrix a 
c	***********************************************************************
	real*8 a(m,m),b(n,n)
	integer n,i,j,k
	if ( m.ne.(n-2) ) stop 'Error in dimensions  (m.ne.n-2) '
	do 1,i=2,n-1
	  do 2,j=2,n-1
	      a(i,j) = b(i,j)  
2	  continue
1	continue 	
	return
	end 
c	***********************************************************************
	subroutine samemsp (a,m,n)
c	store the matrix m in the matrix a 
c	***********************************************************************
	real*4 a(n,n),m(n,n)
	integer n,i,j,k
	do 1,i=2,n-1
	  do 2,j=2,n-1
	      a(i,j) = m(i,j)  
2	  continue
1	continue 	
	do k=1,n
	  a(n,k) = 0.0
	  a(1,k) = 0.0
	  a(k,1) = 0.0
	  a(k,n) = 0.0
	end do
	return
	end 
c	***********************************************************************
	subroutine samevsp (v,w,n)
c	store the vector w in the vector v 
c	***********************************************************************
	real*4 v(n),w(n)
	integer n,i
	do 1,i=2,n-1
	      v(i) = w(i) 
 1	continue
	v(1) = 0.0
	v(n) = 0.0
	return
	end
c	***********************************************************************
	subroutine samev (v,w,n)
c	store the vector w in the vector v 
c	***********************************************************************
	real*8 v(n),w(n)
	integer n,i
	do 1,i=2,n-1
	      v(i) = w(i) 
 1	continue
	v(1) = 0.0d0
	v(n) = 0.0d0
	return
	end
c	***********************************************************************
	subroutine samevcore (v,w,m,n)
c	store the vector w in the vector v 
c	***********************************************************************
	real*8 v(m),w(n)
	integer n,i
	if (m.ne.(n-2)) stop ' Error in dimensions (m.ne.n-2) '
	do 1,i=2,n-1
	      v(i) = w(i) 
 1	continue
	return
	end
!c	***********************************************************************
!	subroutine operaux (a,n,  b,c,d,e, ll,mm,dd,maux1,maux2) 
!c	***********************************************************************
!	real*8 a(n,n),b(n,n),c(n,n),d(n,n),e(n,n) 
!	real*8 maux1(n,n),maux2(n,n),ll(n),mm(n),dd
!	integer n
!	call mulmm(a,c,e,n)
!	call unit(maux1,n)
!	call resmm(maux2,maux1,a,n)			! maux2 = 1 - c e
!	call mynvdpnd(maux2,n,dd,ll,mm)			! maux2 = 1 / (1-ce) 
!	call mulmm(a,c,d,n)
!	call resmm(maux1,b,a,n)				! a = b - c d
!	call mulmm(a,maux2,maux1,n) 			! a = cax2 * (b-cd) 
!	return 
!	end
!c	***********************************************************************
!	subroutine invmcore (a,acore,n,  dd,ll,mm) 
!c	***********************************************************************
!	real*8 a(n,n), acore(n-2,n-2)
!	real*8 ll(n-2),mm(n-2),dd
!	integer i,n,j,k
!
!	do i=2,n-1
!	  do j=2,n-1
!	    acore(i-1,j-1) = a(i,j)
!	  end do
!	end do
!	call mynvdpnd (acore,n-2,dd,ll,mm)		
!	do i=2,n-1
!	  do j=2,n-1
!	    a(i,j) = acore(i-1,j-1)
!	  end do
!	end do
!	do k=1,n
!	  a(1,k) = 0.d0
!	  a(n,k) = 0.d0
!	  a(k,1) = 0.d0
!	  a(k,n) = 0.d0
!	end do
!
!	return 
!	end
c	***********************************************************************
	subroutine mulmv(a,b,c,n)
c	do a(i)=b(i,j)*c(j). a, b, and c must be distint
c	***********************************************************************
	implicit none

	integer n,i,j
	real*8 a(n),b(n,n),c(n),sum

	do 1,i=2,n-1
	  sum=0.0d0
	  do 2,j=2,n-1
	    sum=sum+ (b(i,j)) * (c(j))
 2	  continue
	  a(i)=sum
 1	continue
	a(1) = 0.0d0
	a(n) = 0.0d0
	return
	end



cc	***********************************************************************
	subroutine muymmv(w,b,c,v,n)
c       c(i,j) is diagonall and will be inverted. Let us call Z(i)=c(i,i)^(-1)
c	Z(i) and v(i) are vectors. multiply first Z(i) and v(i) 
c	them multiply b and the previous product. w(i)=b(i,j)*(Z(j)+v(j))
c	***********************************************************************
	real*8 w(n),b(n,n),c(n,n),v(n), sum
	integer n,i,j,k
	do 1,i=2,n-1
	  sum=0.0d0
	  do 2,j=2,n-1
	    sum=sum+ (b(i,j)) * (v(j)/c(j,j))
 2	  continue
	  w(i)=sum
 1	continue
	w(1) = 0.0d0
	w(n) = 0.0d0
	return
	end
cc	***********************************************************************
	subroutine muymvv(w,b,u,v,n)
c       u(i) is to be inverted. Let us call Z=u^(-1)
c	Z(i) and v(i) are vectors. multiply first Z(i) and v(i) 
c	them multiply b and the previous product. w(i)=b(i,j)*(Z(j)+v(j))
c	***********************************************************************
	real*8 w(n),u(n),b(n,n),v(n), sum
	integer n,i,j,k
	do 1,i=2,n-1
	  sum=0.0d0
	  do 2,j=2,n-1
	    sum=sum+ (b(i,j)) * (v(j)/u(j))
 2	  continue
	  w(i)=sum
 1	continue
	w(1) = 0.0d0
	w(n) = 0.0d0
	return
	end
c	***********************************************************************
	subroutine mulmvv(w,b,u,v,n)
c	u(i) and v(i) are vectors. multiply first u(i) and v(i) 
c	them multiply b and the previous product. w(i)=b(i,j)*(u(j)+v(j))
c	***********************************************************************
	real*8 w(n),u(n),b(n,n),v(n), sum
	integer n,i,j,k
	do 1,i=2,n-1
	  sum=0.0d0
	  do 2,j=2,n-1
	    sum=sum+ (b(i,j)) * (u(j)+v(j))
 2	  continue
	  w(i)=sum
 1	continue
	w(1) = 0.0d0
	w(n) = 0.0d0
	return
	end
c	***********************************************************************
	subroutine mulvmv(a,u,b,v,n)
c	u(i) and v(i) are vectors. store u(i) and v(i) in the diagonal 
c	elements of square matrixes. then do a(i,j)= u(i,i)*b(i,j)*v(j,j)
c	***********************************************************************
	real*8 a(n,n),u(n),b(n,n),v(n)
	integer n,i,j,k
	do i=2,n-1
	  do j=2,n-1
	    a(i,j)=(b(i,j)) * (v(j))
	  end do
	end do
	do i=2,n-1
	  do j=2,n-1
	    a(i,j)=(u(i)) * (a(i,j))
	  end do
	end do
	do k=1,n
	  a(1,k) = 0.d0
	  a(n,k) = 0.d0
	  a(k,1) = 0.d0
	  a(k,n) = 0.d0
	end do

	return
	end
c	***********************************************************************
	subroutine mulmm(a,b,c,n)
c	***********************************************************************
	real*8 a(n,n), b(n,n), c(n,n)
	integer n,i,j,k

!	do i=2,n-1
!	  do j=2,n-1
!	    a(i,j)= 0.d00
!	    do k=2,n-1
!		a(i,j) = a(i,j) + b(i,k) * c(k,j)
!	    end do
!	  end do
!	end do
	do j=2,n-1
	   do i=2,n-1
	      a(i,j)=0.d0
	   enddo
	   do k=2,n-1
	      do i=2,n-1
		 a(i,j)=a(i,j)+b(i,k)*c(k,j)
	      enddo
	   enddo
	enddo
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do

	return
	end
c	***********************************************************************
	subroutine summm(a,b,c,n)
c	***********************************************************************
	real*8 a(n,n), b(n,n), c(n,n)
	integer n,i,j,k

	do i=2,n-1
	  do j=2,n-1
	    a(i,j)= b(i,j) + c(i,j)
	  end do
	end do
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do

	return
	end
c	***********************************************************************
	subroutine resmm(a,b,c,n)
c	***********************************************************************
	real*8 a(n,n), b(n,n), c(n,n)
	integer n,i,j,k

	do i=2,n-1
	  do j=2,n-1
	    a(i,j)= b(i,j) - c(i,j)
	  end do
	end do
	do k=1,n
	  a(n,k) = 0.0d0
	  a(1,k) = 0.0d0
	  a(k,1) = 0.0d0
	  a(k,n) = 0.0d0
	end do

	return
	end
c	***********************************************************************
	subroutine mulvv(a,b,c,n)
c	a(i)=b(i)*c(i)
c	***********************************************************************
	real*8 a(n),b(n),c(n)
	integer n,i
	do 1,i=2,n-1
	  a(i)= (b(i)) * (c(i))
 1	continue
	a(1) = 0.0d0
	a(n) = 0.0d0
	return
	end
c	***********************************************************************
	subroutine sumvv(a,b,c,n)
c	a(i)=b(i)+c(i)
c	***********************************************************************
	implicit none

	integer n,i
	real*8 a(n),b(n),c(n)

	do 1,i=2,n-1
	  a(i)= (b(i)) + (c(i))
 1	continue
	a(1) = 0.0d0
	a(n) = 0.0d0
	return
	end

c	***********************************************************************
	subroutine sumvvv(a,b,c,d,n)
c	a(i)=b(i)+c(i)+d(i)
c	***********************************************************************
	real*8 a(n),b(n),c(n),d(n)
	integer n,i
	do 1,i=2,n-1
	  a(i)= b(i) + c(i) + d(i)
 1	continue
	a(1) = 0.0d0
	a(n) = 0.0d0
	return
	end
c	***********************************************************************
	subroutine resvv(a,b,c,n)
c	a(i)=b(i)-c(i)
c	***********************************************************************
	real*8 a(n),b(n),c(n)
	integer n,i
	do 1,i=2,n-1
	  a(i)= (b(i)) - (c(i))
 1	continue
	a(1) = 0.0d0
	a(n) = 0.0d0
	return
	end
c	***********************************************************************
	subroutine zerom(a,n)
c	a(i,j)= 0.0
c	***********************************************************************

	implicit none

	integer n,i,j
	real*8 a(n,n)

	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0d0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zeromsp (a,n)
c	a(i,j)= 0.0
c	***********************************************************************
	real*4 a(n,n)
	integer n,i,j
	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero4m(a,b,c,d,n)
c	a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0 
c	***********************************************************************
	real*8 a(n,n), b(n,n), c(n,n), d(n,n)
	integer n,i,j
	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0d0
	    b(i,j) = 0.0d0
	    c(i,j) = 0.0d0
	    d(i,j) = 0.0d0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero4msp(a,b,c,d,n)
c	a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0 
c	***********************************************************************
	real*4 a(n,n), b(n,n), c(n,n), d(n,n)
	integer n,i,j
	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0
	    b(i,j) = 0.0
	    c(i,j) = 0.0
	    d(i,j) = 0.0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero3m(a,b,c,n)
c	a(i,j) = b(i,j) = c(i,j) = 0.0 
c	**********************************************************************
	real*8 a(n,n), b(n,n), c(n,n)
	integer n,i,j
	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0d0
	    b(i,j) = 0.0d0
	    c(i,j) = 0.0d0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero3msp(a,b,c,n)
c	a(i,j) = b(i,j) = c(i,j) = 0.0 
c	***********************************************************************
	real*4 a(n,n), b(n,n), c(n,n)
	integer n,i,j
	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0
	    b(i,j) = 0.0
	    c(i,j) = 0.0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero2m(a,b,n)
c	a(i,j) = b(i,j) = 0.0 
c	***********************************************************************
	real*8 a(n,n), b(n,n)
	integer n,i,j
	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0d0
	    b(i,j) = 0.0d0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero2msp(a,b,n)
c	a(i,j) = b(i,j) = 0.0 
c	***********************************************************************
	real*4 a(n,n), b(n,n)
	integer n,i,j
	do 1,i=1,n
	  do 2,j=1,n
	    a(i,j) = 0.0
	    b(i,j) = 0.0
 2	  continue
 1	continue
	return
	end
c	***********************************************************************
	subroutine zerov(a,n)
c	a(i)= 0.0
c	***********************************************************************
	real*8 a(n)
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0d0
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero4v(a,b,c,d,n)
c	a(i) = b(i) = c(i) = d(i,j) = 0.0
c	***********************************************************************
	real*8 a(n), b(n), c(n), d(n)
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0d0
	    b(i) = 0.0d0
	    c(i) = 0.0d0
	    d(i) = 0.0d0
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero3v(a,b,c,n)
c	a(i) = b(i) = c(i) = 0.0
c	***********************************************************************
	real*8 a(n), b(n), c(n) 
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0d0
	    b(i) = 0.0d0
	    c(i) = 0.0d0
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero2v(a,b,n)
c	a(i) = b(i) = 0.0
c	***********************************************************************
	real*8 a(n), b(n) 
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0d0
	    b(i) = 0.0d0
 1	continue
	return
	end
c	***********************************************************************
	subroutine zerovsp(a,n)
c	a(i)= 0.0
c	***********************************************************************
	real*4 a(n)
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero4vsp(a,b,c,d,n)
c	a(i) = b(i) = c(i) = d(i) = 0.0
c	***********************************************************************
	real*4 a(n), b(n), c(n), d(n)
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0
	    b(i) = 0.0
	    c(i) = 0.0
	    d(i) = 0.0
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero3vsp(a,b,c,n)
c	a(i) = b(i) = c(i) = 0.0
c	**********************************************************************
	real*4 a(n), b(n), c(n)
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0
	    b(i) = 0.0
	    c(i) = 0.0
 1	continue
	return
	end
c	***********************************************************************
	subroutine zero2vsp(a,b,n)
c	a(i) = b(i) = 0.0
c	***********************************************************************
	real*4 a(n), b(n) 
	integer n,i
	do 1,i=1,n
	    a(i) = 0.0
	    b(i) = 0.0
 1	continue
	return
	end
c	***********************************************************************
	!subroutine zerojt(a,n1,n2,n3)
	subroutine zerovdim3(a,n1,n2,n3)  ! sustituye a zerojt de cristina
c	a(i,j,k)= 0.0 
c	jt(icol,nisos,nb+1,n)
c	***********************************************************************
!	real*4 a(9,34,n)
!	integer n,i,j,k,icol,ic
	real*4 a(n1,n2,n3)
	integer n1,n2,n3,i,j,k

	 do 2,i=1,n1
	  do 3,j=1,n2
	    do 4,k=1,n3
		a(i,j,k) = 0.0
 4	    continue
 3	  continue
 2	 continue

	return
	end
c	***********************************************************************