Commit bea7b78a authored by Serhiy Mochalskyy's avatar Serhiy Mochalskyy
Browse files

coils version

parent e0139aa0
Pipeline #5038 skipped
! ----------------------------------------------------------------------
subroutine bfield_c_par(xp,yp,zp,bx,by,bz,n,x,y,z,phi,ntri)
! ----------------------------------------------------------------------
! 27.06.01
! purpose:
!
!
! ----------------------------------------------------------------------
implicit none
real,dimension (n) :: xp,yp,zp,bx,by,bz
real,dimension (ntri,3) :: x,y,z,phi
integer :: n,ntri
! ----------------------------------------------------------------------
! local variables
real,dimension (n):: x1,y1,z1,x2,y2,z2,x3,y3,z3,sn
real :: s1,s2,s3 &
,d221,d232,d213,al1,al2,al3 &
,vx32,vy32,vz32,vx13,vy13,vz13 &
,ata1,ata2,ata3,v21,v32,v13,at &
,s21,s22,s23,dp1,dm1,dp2,dm2,dp3,dm3 &
,ap1,am1,ap2,am2,ap3,am3 &
,h,ln1,ln2,ln3,ar1,ar2,ar3 &
,x21,y21,z21,x32,y32,z32,x13,y13,z13,vx,vy,vz&
,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3 &
,nx,ny,nz,pi41,area,d21,d32,d13,jx,jy,jz &
,dep1,dep2,dep3,dem1,dem2,dem3
integer :: i,k
! ----------------------------------------------------------------------
pi41 = .125/asin(1.)
bx=0.
by=0.
bz=0.
do i=1,n
do k=1,ntri
x21 = x(k,2) - x(k,1)
y21 = y(k,2) - y(k,1)
z21 = z(k,2) - z(k,1)
x32 = x(k,3) - x(k,2)
y32 = y(k,3) - y(k,2)
z32 = z(k,3) - z(k,2)
x13 = x(k,1) - x(k,3)
y13 = y(k,1) - y(k,3)
z13 = z(k,1) - z(k,3)
d221 = x21**2+y21**2+z21**2
d232 = x32**2+y32**2+z32**2
d213 = x13**2+y13**2+z13**2
d21 = sqrt(d221)
d32 = sqrt(d232)
d13 = sqrt(d213)
nx = -y21*z13 + z21*y13
ny = -z21*x13 + x21*z13
nz = -x21*y13 + y21*x13
area = 1./sqrt(nx*nx+ny*ny+nz*nz)
jx = (x32*phi(k,1)+x13*phi(k,2)+x21*phi(k,3))*area*pi41
jy = (y32*phi(k,1)+y13*phi(k,2)+y21*phi(k,3))*area*pi41
jz = (z32*phi(k,1)+z13*phi(k,2)+z21*phi(k,3))*area*pi41
nx = nx*area
ny = ny*area
nz = nz*area
tx1 = (y32*nz-z32*ny)
ty1 = (z32*nx-x32*nz)
tz1 = (x32*ny-y32*nx)
tx2 = (y13*nz-z13*ny)
ty2 = (z13*nx-x13*nz)
tz2 = (x13*ny-y13*nx)
tx3 = (y21*nz-z21*ny)
ty3 = (z21*nx-x21*nz)
tz3 = (x21*ny-y21*nx)
! ----------------------------------------------------------------------
x1(i) = x(k,1) - xp(i)
y1(i) = y(k,1) - yp(i)
z1(i) = z(k,1) - zp(i)
x2(i) = x(k,2) - xp(i)
y2(i) = y(k,2) - yp(i)
z2(i) = z(k,2) - zp(i)
x3(i) = x(k,3) - xp(i)
y3(i) = y(k,3) - yp(i)
z3(i) = z(k,3) - zp(i)
sn(i) = nx*x1(i)+ny*y1(i)+nz*z1(i)
h = abs(sn(i))
s21 = x1(i)**2+y1(i)**2+z1(i)**2
s22 = x2(i)**2+y2(i)**2+z2(i)**2
s23 = x3(i)**2+y3(i)**2+z3(i)**2
s1 = sqrt(s21)
s2 = sqrt(s22)
s3 = sqrt(s23)
al1 = alog((s2+s1+d21)/(s1+s2-d21))
al2 = alog((s3+s2+d32)/(s3+s2-d32))
al3 = alog((s1+s3+d13)/(s1+s3-d13))
!---------------------------------------------------------------------
ar1 = x1(i)*tx3+y1(i)*ty3+z1(i)*tz3
ar2 = x2(i)*tx1+y2(i)*ty1+z2(i)*tz1
ar3 = x3(i)*tx2+y3(i)*ty2+z3(i)*tz2
dp1 = .5*(s22-s21+d221)
dp2 = .5*(s23-s22+d232)
dp3 = .5*(s21-s23+d213)
dm1 = dp1-d221
dm2 = dp2-d232
dm3 = dp3-d213
ap1 = ar1*dp1
dep1 = ar1**2+h*d221*(h+s2)
ap2 = ar2*dp2
dep2 = ar2**2+h*d232*(h+s3)
ap3 = ar3*dp3
dep3 = ar3**2+h*d213*(h+s1)
am1 = ar1*dm1
dem1 = ar1**2+h*d221*(h+s1)
am2 = ar2*dm2
dem2 = ar2**2+h*d232*(h+s2)
am3 = ar3*dm3
dem3 = ar3**2+h*d213*(h+s3)
ata1 = atan2(ap1*dem1-am1*dep1,dep1*dem1+ap1*am1)
ata2 = atan2(ap2*dem2-am2*dep2,dep2*dem2+ap2*am2)
ata3 = atan2(ap3*dem3-am3*dep3,dep3*dem3+ap3*am3)
at = sign(1.,sn(i))*(ata1+ata2+ata3)
vx = -nx*at + al1*tx3/d21+al2*tx1/d32+al3*tx2/d13
vy = -ny*at + al1*ty3/d21+al2*ty1/d32+al3*ty2/d13
vz = -nz*at + al1*tz3/d21+al2*tz1/d32+al3*tz2/d13
bx(i) = bx(i)+ vy*jz-vz*jy
by(i) = by(i)+ vz*jx-vx*jz
bz(i) = bz(i)+ vx*jy-vy*jx
enddo
enddo
end subroutine bfield_c_par
! ----------------------------------------------------------------------
subroutine bfield_c(xp,yp,zp,bx,by,bz,n,x,y,z,phi,ntri)
! ----------------------------------------------------------------------
......@@ -184,21 +30,6 @@
bx=0.
by=0.
bz=0.
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE( &
!$OMP x1,y1,z1,x2,y2,z2,x3,y3,z3,sn &
!$OMP , s1,s2,s3 &
!$OMP ,d221,d232,d213,al1,al2,al3 &
!$OMP ,vx32,vy32,vz32,vx13,vy13,vz13 &
!$OMP ,ata1,ata2,ata3,v21,v32,v13,at &
!$OMP ,s21,s22,s23,dp1,dm1,dp2,dm2,dp3,dm3 &
!$OMP ,ap1,am1,ap2,am2,ap3,am3 &
!$OMP ,h,ln1,ln2,ln3,ar1,ar2,ar3 &
!$OMP ,x21,y21,z21,x32,y32,z32,x13,y13,z13,vx,vy,vz&
!$OMP ,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3 &
!$OMP ,nx,ny,nz,area,d21,d32,d13,jx,jy,jz &
!$OMP ,dep1,dep2,dep3,dem1,dem2,dem3 &
!$OMP ,i,k )
!$OMP DO
do i=1,n
do k=1,ntri
......@@ -293,7 +124,5 @@
bz(i) = bz(i)+ vx*jy-vy*jx
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end subroutine bfield_c
! ----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine bfield_par_par(b_par, xp, yp, zp, xu, yu, zu, n, n_b, n_e, x, y, z, ntri, j_pot)
! ----------------------------------------------------------------------
! 23.09.09
......@@ -140,11 +140,6 @@ do k=1,ntri
vy = -ny*at + al1*ty3/d21+al2*ty1/d32+al3*ty2/d13
vz = -nz*at + al1*tz3/d21+al2*tz1/d32+al3*tz2/d13
!if(rank==1 .AND. j1==640 ) write(500+rank,*) n_b,n_e
!if(rank==1 .AND. j2==640 ) write(510+rank,*) n_b,n_e
!if(rank==1 .AND. j3==640 ) write(520+rank,*) n_b,n_e
if( (j1>=n_b .AND. j1<=n_e) ) then
......@@ -152,7 +147,6 @@ do k=1,ntri
+ (vz*jx1-vx*jz1)*yu(i) &
+ (vx*jy1-vy*jx1)*zu(i)
b_par(i,j1) = b_par(i,j1)+ temp1
!if(rank==1 .AND. j1==640 ) write(130+rank,*) i, b_par(i,j1), temp1
endif
......@@ -163,8 +157,6 @@ do k=1,ntri
+ (vz*jx2-vx*jz2)*yu(i) &
+ (vx*jy2-vy*jx2)*zu(i)
b_par(i,j2) = b_par(i,j2)+ temp2
!if(rank==1 .AND. j2==640 ) write(140+rank,*) i, b_par(i,j2), temp2
endif
......@@ -175,27 +167,14 @@ do k=1,ntri
+ (vz*jx3-vx*jz3)*yu(i) &
+ (vx*jy3-vy*jx3)*zu(i)
b_par(i,j3) = b_par(i,j3)+ temp3
!if(rank==1 .AND. j3==640 ) write(150+rank,*) i, b_par(i,j3), temp3
endif
! if( (rank==1 .AND. j1==320 .AND. i==1) .OR. &
! (rank==1 .AND. j2==320 .AND. i==1) .OR. &
! (rank==1 .AND. j3==320 .AND. i==1)) write(850,*) b_par(1,320)
! if( (rank==0 .AND. j1==320 .AND. i==1) .OR. &
! (rank==0 .AND. j2==320 .AND. i==1) .OR. &
! (rank==0 .AND. j3==320 .AND. i==1)) write(860,*) b_par(1,320)
enddo
endif
enddo
! if(rank==1 ) write(360+rank,*) b_par(:,640)
end subroutine bfield_par_par
......
......@@ -16,13 +16,13 @@
+(zc(nc,1)-zc(1,1))**2+(xc(nc,2)-xc(1,2))**2 &
+(yc(nc,2)-yc(1,2))**2+(zc(nc,2)-zc(1,2))**2
if(dist.gt.1.e-6) then
write(6,*) dist, xc(nc,1)-xc(1,1),yc(nc,1)-yc(1,1) &
,zc(nc,1)-zc(1,1)
write(6,*)
write(6,*) dist, xc(nc,2)-xc(1,2),yc(nc,2)-yc(1,2) &
,zc(nc,2)-zc(1,2)
write(6,*) 'coil_2d is not closed: stop'
stop
write(6,*) dist, xc(nc,1)-xc(1,1),yc(nc,1)-yc(1,1) &
,zc(nc,1)-zc(1,1)
write(6,*)
write(6,*) dist, xc(nc,2)-xc(1,2),yc(nc,2)-yc(1,2) &
,zc(nc,2)-zc(1,2)
write(6,*) 'coil_2d is not closed: stop'
stop
endif
j=0
do k = 1, nc-1
......
......@@ -109,15 +109,22 @@ $(OBJ_DIR)/matrix_rw.o: \
$(OBJ_DIR)/mod_mpi_v.o $(OBJ_DIR)/mod_sca.o
$(OBJ_DIR)/matrix_cc.o: \
$(OBJ_DIR)/mod_coil2d.o $(OBJ_DIR)/mod_solv.o \
$(OBJ_DIR)/mod_icontr.o
$(OBJ_DIR)/mod_icontr.o $(OBJ_DIR)/mod_tri.o \
$(OBJ_DIR)/mod_tri_w.o $(OBJ_DIR)/mod_mpi_v.o\
$(OBJ_DIR)/mod_sca.o
$(OBJ_DIR)/matrix_cp.o: \
$(OBJ_DIR)/mod_coil2d.o $(OBJ_DIR)/mod_solv.o\
$(OBJ_DIR)/mod_tri_p.o $(OBJ_DIR)/mod_icontr.o
$(OBJ_DIR)/mod_tri_p.o $(OBJ_DIR)/mod_icontr.o\
$(OBJ_DIR)/mod_tri.o $(OBJ_DIR)/mod_tri_b.o \
$(OBJ_DIR)/mod_mpi_v.o
$(OBJ_DIR)/matrix_wc.o: \
$(OBJ_DIR)/mod_coil2d.o $(OBJ_DIR)/mod_solv.o\
$(OBJ_DIR)/mod_tri_w.o $(OBJ_DIR)/mod_icontr.o
$(OBJ_DIR)/mod_tri_w.o $(OBJ_DIR)/mod_icontr.o\
$(OBJ_DIR)/mod_tri.o $(OBJ_DIR)/mod_tri_b.o \
$(OBJ_DIR)/mod_mpi_v.o
$(OBJ_DIR)/matrix_rc.o: \
$(OBJ_DIR)/mod_coil2d.o $(OBJ_DIR)/mod_solv.o
$(OBJ_DIR)/mod_coil2d.o $(OBJ_DIR)/mod_solv.o\
$(OBJ_DIR)/mod_mpi_v.o
$(OBJ_DIR)/matrix_pe.o: \
$(OBJ_DIR)/mod_contr_su.o $(OBJ_DIR)/mod_solv.o\
$(OBJ_DIR)/mod_icontr.o $(OBJ_DIR)/mod_tri_p.o\
......@@ -130,7 +137,7 @@ $(OBJ_DIR)/matrix_ep.o: \
$(OBJ_DIR)/matrix_ec.o: \
$(OBJ_DIR)/mod_contr_su.o $(OBJ_DIR)/mod_solv.o\
$(OBJ_DIR)/mod_icontr.o $(OBJ_DIR)/mod_tri_w.o\
$(OBJ_DIR)/mod_coil2d.o
$(OBJ_DIR)/mod_coil2d.o $(OBJ_DIR)/mod_mpi_v.o
$(OBJ_DIR)/matrix_ew.o: \
$(OBJ_DIR)/mod_contr_su.o $(OBJ_DIR)/mod_solv.o\
$(OBJ_DIR)/mod_icontr.o $(OBJ_DIR)/mod_tri_w.o\
......
......@@ -127,11 +127,11 @@ subroutine input
stop 1
end if
if ( nu_coil >0 ) then
write(outp,*) 'ERROR: nu_coil does not include in this code version. &
The value should be 0'
stop 1
end if
! if ( nu_coil >0 ) then
! write(outp,*) 'ERROR: nu_coil does not include in this code version. &
! The value should be 0'
! stop 1
! end if
......
......@@ -38,7 +38,7 @@ program starwall
! B_perp on wall is zero (ideal conductor)
! Check if the amount of MPI tasks more than array dimencions
call control_array_distribution
call control_array_distribution
! Main solver
call solver
......
......@@ -8,30 +8,55 @@
use icontr
use coil2d
use solv
use tri
use tri_w
use mpi_v
use sca
! ----------------------------------------------------------------------
implicit none
integer :: i,i0,i1,ip,ip0,ip1,j,jp,k,ier
integer :: i,i0,i1,ip,ip0,ip1,j,jp,k
real :: pi,mu0
real :: dx1,dy1,dz1,dx2,dy2,dz2,dx3,dy3,dz3
real ,allocatable :: jx(:),jy(:),jz(:)
real ,allocatable :: dima(:,:)
real :: dima_sca,dima_sca2,dima_sum,num
integer j_loc,i_loc,il,jl
logical inside_j,inside_i
! ----------------------------------------------------------------------
print *, 'compute matrix_cc'
if (rank==0) write(6,*) 'compute matrix_cc'
! ----------------------------------------------------------------------
write(6,*) 'ma_cc ntri_c=',ntri_c,'ncoil=',ncoil
allocate(dima(ntri_c,ntri_c),jx(ntri_c),jy(ntri_c),jz(ntri_c),stat=ier)
dima=0.; jx=0.; jy=0.; jz=0.
if (rank==0) write(6,*) 'ma_cc ntri_c=',ntri_c,'ncoil=',ncoil
call tri_induct(x_coil,y_coil,z_coil,ntri_c &
,x_coil,y_coil,z_coil,ntri_c,dima)
allocate(jx(ntri_c),jy(ntri_c),jz(ntri_c),stat=ier)
if (ier /= 0) then
print *,'matrix_cc.f90: Can not allocate local arrays jx, ..., jz'
stop
endif
jx=0.; jy=0.; jz=0.
allocate(nx(ntri_c),ny(ntri_c),nz(ntri_c),tx1(ntri_c), &
ty1(ntri_c),tz1(ntri_c),tx2(ntri_c),ty2(ntri_c), &
tz2(ntri_c),tx3(ntri_c),ty3(ntri_c),tz3(ntri_c), &
d221(ntri_c),d232(ntri_c),d213(ntri_c),area(ntri_c), &
xm(7*ntri_c), ym(7*ntri_c), zm(7*ntri_c),stat=ier)
if (ier /= 0) then
print *,'matrix_cc.f90: Can not allocate local arrays nx, ..., zm'
stop
endif
nx=0.; ny=0.; nz=0.; tx1=0.; ty1=0.;tz1=0.; tx2=0.
ty2=0.; tz2=0.; tx3=0.; ty3=0.; tz3=0.; d221=0.
d232=0.; d213=0.;area=0.; xm=0.; ym=0.; zm=0.
!-----------------------------------------------------------------------
pi = 2.*asin(1.)
mu0 = 4.e-07*pi
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(j)
!$OMP DO
pi = 2.*asin(1.)
mu0 = 4.e-07*pi
do j=1,ntri_c
dx1 = x_coil(j,2)-x_coil(j,3)
dy1 = y_coil(j,2)-y_coil(j,3)
......@@ -42,42 +67,88 @@
dx3 = x_coil(j,1)-x_coil(j,2)
dy3 = y_coil(j,1)-y_coil(j,2)
dz3 = z_coil(j,1)-z_coil(j,2)
jx(j) = dx1*phi_coil(j,1)+dx2*phi_coil(j,2)+dx3*phi_coil(j,3)
jy(j) = dy1*phi_coil(j,1)+dy2*phi_coil(j,2)+dy3*phi_coil(j,3)
jz(j) = dz1*phi_coil(j,1)+dz2*phi_coil(j,2)+dz3*phi_coil(j,3)
enddo
!$OMP END DO
!$OMP END PARALLEL
i0 = 0
jx(j) = dx1*phi_coil(j,1)+dx2*phi_coil(j,2)+dx3*phi_coil(j,3)
jy(j) = dy1*phi_coil(j,1)+dy2*phi_coil(j,2)+dy3*phi_coil(j,3)
jz(j) = dz1*phi_coil(j,1)+dz2*phi_coil(j,2)+dz3*phi_coil(j,3)
call tri_induct_1(j,x_coil,y_coil,z_coil,ntri_c)
call tri_induct_2(ntri_c,j,x_coil,y_coil,z_coil)
enddo
i0 = 0
do j =1,ncoil
do i =1,jtri_c(j)
i1 = i0 +i
ip0 = 0
do jp = 1,ncoil
do ip = 1,jtri_c(jp)
ip1 = ip0 +ip
a_ww(jp,j) = a_ww(jp,j)+.5*(dima(ip1,i1)+dima(i1,ip1)) &
*(jx(ip1)*jx(i1) &
+jy(ip1)*jy(i1) &
+jz(ip1)*jz(i1) )
enddo
ip0 = ip0+jtri_c(jp)
enddo
enddo
call ScaLAPACK_mapping_j(j,j_loc,inside_j)
if (inside_j == .true.) then! If not inside
do i =1,jtri_c(j)
i1 = i0 +i
ip0 = 0
do jp = 1,ncoil
call ScaLAPACK_mapping_i(jp,i_loc,inside_i)
if (inside_i == .true.) then! If not inside
do ip = 1,jtri_c(jp)
ip1 = ip0 +ip
dima_sca=0
dima_sca2=0
call tri_induct_3_2(ntri_c,ntri_c,ip1,i1,x_coil,y_coil,z_coil,dima_sca)
call tri_induct_3_2(ntri_c,ntri_c,i1,ip1,x_coil,y_coil,z_coil,dima_sca2)
dima_sum=dima_sca+dima_sca2
a_ww_loc(i_loc,j_loc) = a_ww_loc(i_loc,j_loc)+&
.5*dima_sum &
*(jx(ip1)*jx(i1) &
+jy(ip1)*jy(i1) &
+jz(ip1)*jz(i1) )
enddo
endif
ip0 = ip0+jtri_c(jp)
enddo
enddo
endif
i0 = i0+jtri_c(j)
enddo
deallocate(jz,jy,jx,dima)
deallocate(jz,jy,jx)
deallocate(nx,ny,nz,tx1,ty1,tz1,tx2,ty2, &
tz2,tx3,ty3,tz3,d221,d232,d213,area,xm,ym,zm)
!-------------------------------------------------------------------
CALL DESCINIT(DESCA, npot_w+ncoil, npot_w+ncoil, NB, NB, 0, 0, CONTEXT, LDA_ww, INFO_A )
if(INFO_A .NE. 0) then
write(6,*) "Something is wrong in matrix_cc CALL DESCINIT DESCA, INFO_A=",INFO_A
stop
endif
if(rank==0) then
open(500,file='coil_inductances.txt',status='unknown')
write(500,*) 'inductances L_ik microHenry'
write(500,*) ' i k L_ik '
endif
do j=1,ncoil
do jp=1,ncoil
write(500,5000) j,jp,mu0*1.e+6*a_ww(j,jp)
5000 format(2i4,1pe12.4)
CALL pdelget('A','D',num, a_ww_loc,j,jp,DESCA)
if(rank==0) write(500,5000) j,jp,mu0*1.e+6*num
5000 format(2i4,1pe12.4)
enddo
enddo
close(500)
if(rank==0) close(500)
!-------------------------------------------------------------------
print *, 'matrix_cc done'
if(rank==0) write(*,*) 'matrix_cc done'
if(rank==0) write(*,*) '==============================================================='
end subroutine matrix_cc
......@@ -9,29 +9,76 @@
use coil2d
use tri_p
use solv
use tri
use tri_b
use mpi_v
! ----------------------------------------------------------------------
implicit none
real ,dimension(:,:),allocatable :: dima,dimb,dxp,dyp,dzp
include "mpif.h"
real ,dimension(:,:),allocatable :: dxp,dyp,dzp
real ,dimension(: ),allocatable :: jx,jy,jz,jx_p,jy_p,jz_p
real :: dx1,dy1,dz1,dx2,dy2,dz2
real :: dx3,dy3,dz3,temp
integer :: i,ic,ic0,ic1,ip,j,k,ier
integer :: i,ic,ic0,ic1,ip,j,k,jp
real :: dima_sum=0,dima_sca=0, dima_sca_b=0
integer :: j_loc,i_loc
logical :: inside_j,inside_i
! ----------------------------------------------------------------------
write(*,*) 'compute matrix_cp'
if(rank==0) write(6,*) 'compute matrix_cp ', 'ntri_p=', ntri_p, 'ntri_c=', ntri_c
! ----------------------------------------------------------------------
allocate(dima(ntri_p,ntri_c),dimb(ntri_c,ntri_p) &
,dxp(ntri_p,3),dyp(ntri_p,3),dzp(ntri_p,3) &
allocate(dxp(ntri_p,3),dyp(ntri_p,3),dzp(ntri_p,3) &
,jx(ntri_c), jy(ntri_c), jz(ntri_c) &
,jx_p(ntri_p),jy_p(ntri_p),jz_p(ntri_p),stat=ier)
if (ier /= 0) then
print *,'Matrix_cp.f90: Can not allocate local arrays nx, ..., zm'
stop
endif
dima=0.; dimb=0.; jx=0.; jy=0.; jz=0.; jx_p=0.; jy_p=0.; jz_p=0.
jx=0.; jy=0.; jz=0.; jx_p=0.; jy_p=0.; jz_p=0.
!-----------------------------------------------------------------------
! ----------------------------------------------------------------------
! Additional arrays that will be used for tri_induct subroutine
! For the production run with n_tri_w=5*10^5 the size of all array will be around 300Mb
! if they will not be distributed among MPI tasks
allocate(nx(ntri_c),ny(ntri_c),nz(ntri_c),tx1(ntri_c), &
ty1(ntri_c),tz1(ntri_c),tx2(ntri_c),ty2(ntri_c), &
tz2(ntri_c),tx3(ntri_c),ty3(ntri_c),tz3(ntri_c), &
d221(ntri_c),d232(ntri_c),d213(ntri_c),area(ntri_c), &
xm(7*ntri_p), ym(7*ntri_p), zm(7*ntri_p), stat=ier)
if (ier /= 0) then
print *,'Matrix_cp.f90: Can not allocate local arrays nx, ..., zm'
stop
endif
nx=0.; ny=0.; nz=0.; tx1=0.; ty1=0.;tz1=0.; tx2=0.
ty2=0.; tz2=0.; tx3=0.; ty3=0.; tz3=0.; d221=0.