39 SUBROUTINE ini_fxbody(FXBIPM, FXBRPM, FXBNOD, FXBGLM,FXBCPM,
40 . FXBCPS, FXBLM, FXBFLS, FXBDLS,FXBMOD,
41 . ITAB, X ,MS, IN, FXB_MATRIX,
42 . FXB_MATRIX_ADD,FXB_LAST_ADRESS,ICODE,NOM_OPT)
49 USE format_mod ,
ONLY : fmw_10i
53#include "implicit_f.inc"
67 INTEGER ITAB(*),FXBIPM(NBIPM,*),FXBNOD(*),FXB_MATRIX_ADD(4,*),FXB_LAST_ADRESS(*),ICODE(*),NOM_OPT(LNOPT1,*)
68 my_real X(3,*),MS(*),IN(*),
69 . FXBRPM(*), FXBGLM(*), FXBCPM(*), FXBCPS(*),
70 . fxblm(*), fxbfls(*), fxbdls(*), fxbmod(*),
75 INTEGER NFX,ID,IDMAST,NMOD,NMST,NBNO,,NTR,ADRGLM,
76 . ADRCP,ADRLM,ADRFLS,ADRDLS,ADRVAR,,ADRMCD,ADRMOD,IMOD,INO,I,LEN,
77 . NLIG,NRES,ILIG,ADRCP2,IR,ADRNOD,NUMNO(10),ISHELL,
78 . ic,j,info,ianim, imin, imax, ref, err, k, l, idn, i1,i2, j2, ref2,j1,size_max,idum1,idum2,idum3,
79 . idum4,flag,ic1,ic2,bcs(6),size_stiff,adr_stiff,il1,il2,idof1,idof2,nmod_max,nmst_f,nmr,
80 . adr_mass,size_mass,reloop,adrmod0,adrlm0,adrfls0,lencp0,lenlm0,lenfls0,lendls0,lenvar0,
81 . nsav_modes,nsav_modesn,nn
83 . freq,beta,omega,dtc1,dtc2,bid,xm(3),zz,rho_invmk,unsn,
norm,fac,fac2,mstot,tole,
84 . kt,kr,rr,skew(9),rdum1,min_diag_stif,fac2x,fac2y,fac2z,dt_mini,dt_cst,omega_min,
85 . inmin,inmax,omega_max,max_freq,
alpha,freq_range,min_cut_off_freq
86 CHARACTER(LEN=NCHARTITLE) :: TITR
87 CHARACTER :: MESS*40,MESS1*40,NWLINE*100,FXBFILE*100
88 INTEGER,
DIMENSION(:,:),
ALLOCATABLE :: TABSL,ITAG_DOF
89 INTEGER,
DIMENSION(:),
ALLOCATABLE :: INDEX,LISTN
90 INTEGER IFLAGI1,IFLAGDBL,IRB,CPT,REST,NLIGNE,M,INFO2,LWORK,NBDYN,NDDL,IM,NMOD0,NMST0
92 INTEGER RPM_L,ADRMODE_L,ADRLM_L,ADRFLS_L,ADRGLM_L,ADRCP_L,ADRCP2_L
94 INTEGER NBNOD,IROT,IDAMP,IBLO,IFILE,,PRINTOUT
95 my_real,
DIMENSION(:),
ALLOCATABLE :: vt,vbid
96 my_real,
DIMENSION(:,:),
ALLOCATABLE :: mass,stif,modes,modest,modes_rb,modes_rbt,stift,
97 . mtemp,mtemp2,mtemp3,tt,t,vectr,modes_sav
99 DOUBLE PRECISION WORK1(1000)
100 DOUBLE PRECISION,
DIMENSION(:),
ALLOCATABLE :: WORK,WR,WI
101 DOUBLE PRECISION,
DIMENSION(:,:),
ALLOCATABLE :: EIG_R,EIG_L,INVMK
103 INTEGER :: LEN_TMP_NAME
104 CHARACTER(len=2148) :: TMP_NAME
112 DATA mess/
'FLEXIBLE BODY : NODES '/
113 DATA mess1/
'FLEXIBLE BODY DEFINITION '/
117 adrmod = fxb_last_adress(1)
118 adrglm = fxb_last_adress(2)
119 adrcp = fxb_last_adress(3)
120 adrcp2 = fxb_last_adress(3)
121 adrlm = fxb_last_adress(4)
122 adrfls = fxb_last_adress(5)
123 adrdls = fxb_last_adress(6)
124 adrvar = fxb_last_adress(7)
125 adrrpm = fxb_last_adress(8)
126 adrmcd = fxb_last_adress(9)
130 min_cut_off_freq = 1000.0 / fac_time
134 adrnod = fxbipm(6,nfx)
136 IF (fxbipm(41,nfx)==2)
THEN
138 IF (printout == 0)
WRITE(iout,1000)
140 CALL fretitl2(titr,nom_opt(lnopt1-ltitr+1,nfx),ltitr)
143 nbnod = fxbipm(3,nfx)
146 adrnod = fxbipm(6,nfx)
147 ishell = fxbipm(16,nfx)
149 iblo = fxbipm(28,nfx)
150 ifile = fxbipm(29,nfx)
151 ianim = fxbipm(36,nfx)
152 size_stiff = fxbipm(42,nfx)
153 size_mass = fxbipm(43,nfx)
154 adr_stiff = fxbipm(44,nfx)
155 adr_mass = fxbipm(45,nfx)
157 size_max =
max(nmod,nme)
159 ALLOCATE(itag_dof(6,nbnod))
160 ALLOCATE(stif(nddl,nddl),mass(nddl,nddl))
161 ALLOCATE(modes_rb(nddl,nme),modes_rbt(nme,nddl))
162 ALLOCATE(modest(nmod,nddl),modes(nddl,nmod))
163 ALLOCATE(vectr(nddl,6),mtemp(nddl,size_max),mtemp2(size_max,size_max))
164 ALLOCATE(mtemp3(size_max,size_max),vt(nddl))
170 fxbipm(14,nfx) = adrrpm
177 fxbrpm(adrrpm+i) = skew(i)
184 IF (fxbnod(adrnod+i-1) < 0)
THEN
185 fxbnod(adrnod+i-1) = abs(fxbnod(adrnod+i-1))
188 fxbipm(18,nfx) = nbnod
195 stif(1:nddl,1:nddl) = zero
197 i1 = fxb_matrix_add(1,adr_stiff+i-1)
198 i2 = fxb_matrix_add(2,adr_stiff+i-1)
199 idof1 = fxb_matrix_add(3,adr_stiff+i-1)
200 idof2 = fxb_matrix_add(4,adr_stiff+i-1)
201 stif(6*(i1-1)+idof1,6*(i2-1)+idof2) = fxb_matrix(adr_stiff+i-1)
202 stif(6*(i2-1)+idof2,6*(i1-1)+idof1) = fxb_matrix(adr_stiff+i-1)
203 IF ((6*(i2-1)+idof2)==(6*(i1-1)+idof1)) min_diag_stif =
min(min_diag_stif,fxb_matrix(adr_stiff+i-1))
210 mass(1:nddl,1:nddl) = zero
212 IF (size_mass == 0)
THEN
215 idn = fxbnod(adrnod+i-1)
216 fxbnod(adrnod+i-1) = idn
217 mass(6*(i-1)+1,6*(i-1)+1) = ms(idn)
218 mass(6*(i-1)+2,6*(i-1)+2) = ms(idn)
219 mass(6*(i-1)+3,6*(i-1)+3) = ms(idn)
220 mass(6*(i-1)+4,6*(i-1)+4) = in(idn)
221 mass(6*(i-1)+5,6*(i-1)+5) = in(idn)
226 i1 = fxb_matrix_add(1,adr_mass+i-1)
227 i2 = fxb_matrix_add(2,adr_mass+i-1)
228 idof1 = fxb_matrix_add(3,adr_mass+i-1)
229 idof2 = fxb_matrix_add(4,adr_mass+i-1)
230 mass(6*(i1-1)+idof1,6*(i2-1)+idof2) = fxb_matrix(adr_mass+i-1)
231 mass(6*(i2-1)+idof2,6*(i1-1)+idof1) = fxb_matrix(adr_mass+i-1)
235 idn = fxbnod(adrnod+i-1)
236 fxbnod(adrnod+i-1) = idn
237 mass(6*(i-1)+1,6*(i-1)+1) = mass(6*(i-1)+1,6*(i-1)+1)+ms(idn)
238 mass(6*(i-1)+2,6*(i-1)+2) = mass(6*(i-1)+2,6*(i-1)+2)+ms(idn)
239 mass(6*(i-1)+3,6*(i-1)+3) = mass(6*(i-1)+3,6*(i-1)+3)+ms(idn)
240 mass(6*(i-1)+4,6*(i-1)+4) = mass(6*(i-1)+4,6*(i-1)+4)+in(idn)
241 mass(6*(i-1)+5,6*(i-1)+5) = mass(6*(i-1)+5,6*(i-1)+5)+in(idn)
242 mass(6*(i-1)+6,6*(i-1)+6) = mass(6*(i-1)+6,6*(i-1)+6)+in
250 fxbipm(7,nfx) = adrmod
257 idn = fxbnod(adrnod+i-1)
258 xm(1) = xm(1) + x(1,idn)*ms(idn)
259 xm(2) = xm(2) + x(2,idn)*ms(idn)
260 xm(3) = xm(3) + x(3,idn)*ms(idn)
261 mstot = mstot + ms(idn)
264 xm(i) = xm(i)/
max(em20,mstot)
267 modes_rb(1:nddl,1:nme) = zero
268 modes_rbt(1:nme,1:nddl) = zero
274 idn = fxbnod(adrnod+k-1)
275 modes_rb(6*(k-1)+j,3*(i-1)+j) = x(i,idn) - xm(i)
282 idn = fxbnod(adrnod+k-1)
283 modes_rb(6*(k-1)+i,9+i) = 1-(x(1,idn)-xm(1))-(x(2,idn)-xm(2))-(x(3,idn)-xm(3))
290 modes_rb(6*(k-1)+3+i,12+i) = one
297 modes_rbt(i,j) = modes_rb(j,i)
303 fxbmod(adrmod) = modes_rb(j,i)
309 fac2x=prscal(modes_rb(1,1), modes_rb(1,1), nddl, mass)
310 fac2y=prscal(modes_rb(1,4), modes_rb(1,4), nddl, mass)
311 fac2z=prscal(modes_rb(1,7), modes_rb(1,7), nddl, mass)
312 fac2 =
max(fac2x,fac2y,fac2z)
327 modes(1:nddl,1:nmod) = zero
328 modest(1:nmod,1:nddl) = zero
329 itag_dof(1:6,1:nbnod) = 0
333 il1 = fxb_matrix_add(1,adr_stiff+i-1)
334 il2 = fxb_matrix_add(2,adr_stiff+i-1)
335 idof1 = fxb_matrix_add(3,adr_stiff+i-1)
336 idof2 = fxb_matrix_add(4,adr_stiff+i-1)
338 itag_dof(idof1,il1) = 1
339 itag_dof(idof2,il2) = 1
344 idn = fxbnod(adrnod+i-1)
350 bcs(2) = (ic1-4*bcs(1))/2
351 bcs(3) = ic1-4*bcs(1)-2*bcs(2)
353 bcs(5) = (ic2-4*bcs(4))/2
354 bcs(6) = ic2-4*bcs(4)-2*bcs(5)
358 IF ((bcs(k)==0).AND.(itag_dof(k,i)>0))
THEN
360 modes(6*(i-1)+k,nmst) = one
374 idn = fxbnod(adrnod+i-1)
382 vectr(cpt+j,5)=x(3,idn)-xm(3)
383 vectr(cpt+j,6)=-(x(2,idn)-xm(2))
386 vectr(cpt+j,4)=-(x(3,idn)-xm(3))
387 vectr(cpt+j,6)=x(1,idn)-xm(1)
390 vectr(cpt+j,4)=x(2,idn)-xm(2)
391 vectr(cpt+j,5)=-(x(1,idn)-xm(1))
405 IF (nmst > 0)
CALL orthrg(vectr,mass,nddl,nmr)
409 IF ((iblo==0).AND.(nmst > 0))
THEN
411 CALL orthsr(modes,vectr,mass,nddl,nmst,nmr)
445 modes(1:nddl,1:nmod) = modes_sav(1:nddl,1:nmod)
456 CALL orthst(modes,mass,nddl,nmst,nmst_f,tole)
466 lencp =lencp -ntr*nmod0*nme+ntr*nmod*nme
467 lenlm =lenlm -nmod0+nmod
468 lenfls=lenfls-nmst0*(2*nmod0-nmst0+1)/2+nmst*(2*nmod-nmst+1)/2
469 lendls=lendls-nmod0+nmst0+nmod+nmst
470 lenvar=lenvar-nmod0+nmod
471 fxbipm(38,nfx)=
min(nmod,fxbipm(38,nfx))
475 modest(j,i)=modes(i,j)
481 fxbmod(adrmod) = fac*modes(j,i)
490 fxbipm(10,nfx) = adrlm
492 mtemp(1:nddl,1:nmod) = matmul(mass(1:nddl,1:nddl),modes(1:nddl,1:nmod))
493 mtemp2(1:nmod,1:nmod) = matmul(modest(1:nmod,1:nddl),mtemp(1:nddl,1:nmod))
496 fxblm(adrlm) = fac2*mtemp2(i,i)
504 fxbipm(11,nfx) = adrfls
506 mtemp(1:nddl,1:nmst) = matmul(stif(1:nddl,1:nddl),modes(1:nddl,1:nmst))
507 mtemp3(1:nmst,1:nmst) = matmul(modest(1:nmst,1:nddl),mtemp(1:nddl,1:nmst))
513 fxbfls(adrfls) = fac2*mtemp3(i,j)
523 ALLOCATE(wi(nmst),wr(nmst),eig_r(nmst,nmst),eig_l(nmst,nmst),invmk(nmst,nmst))
528 mtemp2(i,j) = one/
max(em20,mtemp2(i,j))
535 invmk(1:nmst,1:nmst) = matmul(mtemp2(1:nmst,1:nmst),mtemp3(1:nmst,1:nmst))
537#ifndef WITHOUT_LINALG
539 CALL dgeev(
'N',
'V',nmst,invmk,nmst,wr,wi,eig_l,nmst,eig_r,nmst,work1,lwork,info2)
540 lwork =
max( 1000, int( work1( 1 ) ) )
541 ALLOCATE(work(lwork))
543 CALL dgeev(
'N',
'V',nmst,invmk,nmst,wr,wi,eig_l,nmst,eig_r,nmst,work,lwork,info2)
549 WRITE(*,*)
'The algorithm failed to compute eigenvalues.'
556 IF (wi(i)==zero)
THEN
557 omega_min =
min(omega_min,sqrt(wr(i)))
561 omega_max =
max(freq_range*omega_min,two*pi*min_cut_off_freq)
562 dt_mini = two/omega_max
565 IF (wi(i)==zero)
THEN
566 rho_invmk =
max(rho_invmk,abs(wr(i)))
567 IF ((two/sqrt(abs(wr(i))) > dt_mini))
THEN
568 nsav_modesn = nsav_modesn+1
573 IF ((nmst - nsav_modesn > 0).AND.(reloop==0))
THEN
574 nsav_modes = nsav_modesn
575 ALLOCATE(modes_sav(nddl,nsav_modes))
576 modes_sav(1:nddl,1:nsav_modes) = zero
579 IF ((two/sqrt(abs(wr(i))) > dt_mini).AND.((wi(i)==zero)))
THEN
583 modes_sav(j,nn)=modes_sav(j,nn)+modes(j,k)*eig_r(k,i)
590 DEALLOCATE(wi,wr,eig_r,eig_l,work,invmk)
592 fxbrpm(fxbipm(14,nfx))=zep9*two/sqrt(rho_invmk)
594 IF (nmst - nsav_modesn == 0)
EXIT
602 IF (
ALLOCATED(modes_sav))
DEALLOCATE (modes_sav)
612 fxbrpm(adrrpm+1)=zero
616 fxbrpm(adrrpm+1)=zero
623 fxbipm(8,nfx) = adrglm
625 mtemp(1:nddl,1:nme) = matmul(mass(1:nddl,1:nddl),modes_rb(1:nddl,1:nme))
626 mtemp2(1:nme,1:nme) = matmul(modes_rbt(1:nme,1:nddl),mtemp(1:nddl,1:nme))
632 fxbglm(adrglm) = mtemp2(i,j)
642 fxbipm(9,nfx) = adrcp
651 vt(cpt+i) = modes(cpt+j,k)
652 vt(cpt+3+i) = modes(cpt+3+j,k)
656 mtemp(l,k)=fac*prscal(vt,modes_rb(1,l),nddl,mass)
662 fxbcpm(adrcp) = mtemp(l,k)
674 fxbipm(9,nfx) = adrcp2
683 vt(cpt+i) = modes_rb(cpt+j,k)
684 vt(cpt+3+i) = modes_rb(cpt+3+j,k)
688 mtemp(k,l)=fac*prscal(modes(1,l),vt,nddl,stif)
694 fxbcps(adrcp2) = mtemp(l,k)
706 fxbipm(13,nfx) = adrvar
707 adrvar=adrvar+nmod+nme
708 fxbipm(15,nfx) = adrmcd
709 adrmcd=adrmcd+nme*nme
715 WRITE(iout,1100) id,trim(titr),nmst,nme,nbnod,fxbrpm(fxbipm(14,nfx))
717 WRITE(iout,fmt=fmw_10i)(itab(fxbnod(adrnod+i-1)),i=1,nbnod)
723 rpm_l = fxbipm(14,nfx)
724 adrmode_l = fxbipm(7,nfx)
725 adrlm_l = fxbipm(10,nfx)
726 adrfls_l = fxbipm(11,nfx)
727 adrglm_l = fxbipm(8,nfx)
728 adrcp_l = fxbipm(9,nfx)
729 adrcp2_l = fxbipm(9,nfx)
740 OPEN(unit=ref,file=tmp_name(1:len_tmp_name),
741 . access=
'SEQUENTIAL',form=
'FORMATTED',status=
'UNKNOWN')
743 WRITE(ref,
'(A)')
'# FLEXIBLE BODY DATA'
744 WRITE(ref,
'(A)')
'#--1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|--10---|'
745 WRITE(ref,
'(A)')
'# Param'
746 WRITE(ref,
'(A)')
"# Nbmod Nbstat Nbnod Irot Idamp Iblo Ifile Intype"
748 WRITE(ref,
'(7I8)') nmod,nmst,nbnod,ishell,idamp,iblo,ifile
749 WRITE(ref,
'(A)')
'# Nodes'
755 WRITE(ref,
'(10I8)') (itab(fxbnod(adrnod+cpt+k-1)),k=1,10)
758 rest = nbnod-10*nligne
759 IF (rest
WRITE(ref
'(10I8)') (itab(fxbnod(adrnod+cpt+k-1)),k=1,rest)
761 max_freq = half*omega_max/pi
762 WRITE(ref
'(A)')
'# Mrot11 Mrot12 Mrot13 Mrot21 Mrot22'
763 WRITE(ref,
'(5(1PE16.9))') fxbrpm
764 WRITE(ref
'(A)''# Mrot23 Mrot31 Mrot32 Mrot33 Freq'
765 WRITE(ref,
'(5(1PE16.9))') fxbrpm(rpm_l+6),fxbrpm(rpm_l+7),fxbrpm(rpm_l+8),fxbrpm(rpm_l+9),max_freq
770 WRITE(ref,
'(A)')
'# RAYLEIGH DAMPING'
773 WRITE(ref,
'(2(1PE16.9))') fxbrpm(adrrpm-4),fxbrpm(adrrpm-3)
776 WRITE(ref,
'(A)')
'# BLOCK5 - RB MODES'
783 WRITE(ref,
'(A)')
'# X Y Z XX YY'
784 WRITE(ref,
'(A)')
'# ZZ'
786 WRITE(ref,
'(5(1PE16.9))') (fxbmod(adrmode_l+k-1),k=1,5)
787 WRITE(ref,
'((1PE16.9))') fxbmod(adrmode_l+6-1)
788 adrmode_l = adrmode_l + 6
794 WRITE(ref,
'(A)')
'# BLOCK6 - REDUCED MODES ROTATION '
795 WRITE(ref,
'(A)')
'# X Y Z XX YY'
796 WRITE(ref,
'(A)')
'# ZZ'
800 WRITE(ref,
'(A)')
'# BLOCK7 - REDUCED MODES'
803 WRITE(ref,
'(A)')
'# X Y Z XX YY'
804 WRITE(ref,
'(A)')
'# ZZ'
806 WRITE(ref,
'(5(1PE16.9))') (fxbmod(adrmode_l+k-1),k=1,5)
807 WRITE(ref,
'((1PE16.9))') fxbmod(adrmode_l+6-1)
808 adrmode_l = adrmode_l + 6
814 WRITE(ref,
'(A)')
'# BLOCK8 - REDUCED DIAG MASS MATRIX'
815 WRITE(ref,
'(A)')
'# X Y Z XX YY'
816 WRITE(ref,
'(A)')
'# ZZ'
821 WRITE(ref,
'(5(1PE16.9))') (fxblm(adrlm_l+k-1),k=1,5)
822 adrlm_l = adrlm_l + 5
826 WRITE(ref,
'(5(1PE16.9))') (fxblm(adrlm_l+k-1),k=1,rest)
827 adrlm_l = adrlm_l + 1
832 WRITE(ref,
'(A)')
'# BLOCK9 - REDUCED FULL STIFF MATRIX'
833 WRITE(ref,
'(A)')
'# X Y Z XX YY'
834 WRITE(ref,
'(A)')
'# ZZ'
836 size_max = ((nmst+1)*nmst)/2
840 WRITE(ref,
'(5(1PE16.9))') (fxbfls(adrfls_l+k-1),k=1,5)
841 adrfls_l = adrfls_l + 5
843 rest = size_max-5*nligne
845 WRITE(ref,
'(5(1PE16.9))') (fxbfls(adrfls_l+k-1),k=1,rest)
846 adrfls_l = adrfls_l + 1
851 WRITE(ref,
'(A)')
'# BLOCK11 - RB MODES RB MATRIX'
852 WRITE(ref,
'(A)')
'# X Y Z XX YY'
853 WRITE(ref,
'(A)')
'# ZZ'
855 size_max = ((nme+1)*nme)/2
859 WRITE(ref,
'(5(1PE16.9))') (fxbglm(adrglm_l+k-1),k=1,5)
860 adrglm_l = adrglm_l + 5
862 rest = size_max-5*nligne
864 WRITE(ref
'(5(1PE16.9))') (fxbglm(adrglm_l+k-1),k=1,rest)
865 adrglm_l = adrglm_l + 1
870 WRITE(ref,
'(A)')
'# BLOCK12 - COUPLE MASS PROJECTION'
871 WRITE(ref,
'(A)')
'# X Y Z XX YY'
872 WRITE(ref,
'(A)')
'# ZZ'
875 WRITE(ref,
'(A,7I8)')
'#',l
880 WRITE(ref,
'(5(1PE16.9))') (fxbcpm(adrcp_l+k-1),k=1,5)
881 adrcp_l = adrcp_l + 5
883 rest = size_max-5*nligne
885 WRITE(ref,
'(5(1PE16.9))') (fxbcpm(adrcp_l+k
886 adrcp_l = adrcp_l + 1
892 WRITE(ref,
'(A)')
'# BLOCK13 - COUPLE STIFF PROJECTION'
893 WRITE(ref,
'(A)')
'# X Y Z XX YY'
894 WRITE(ref,
'(A)')
'# ZZ'
897 WRITE(ref,
'(A,7I8)')
'#',l
902 WRITE(ref,
'(5(1PE16.9))') (fxbcps(adrcp2_l+k-1),k=1,5)
903 adrcp2_l = adrcp2_l + 5
905 rest = size_max-5*nligne
907 WRITE(ref,
'(5(1PE16.9))') (fxbcps(adrcp2_l+k-1),k=1,rest)
908 adrcp2_l = adrcp2_l + 1
912 WRITE(ref,
'(A)')
'#--1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|--10---|'
914 CLOSE(unit=ref,status=
'KEEP')
922 DEALLOCATE(itag_dof,stif,mass,modes_rb,modes_rbt)
923 DEALLOCATE(modest,modes,vectr,mtemp,mtemp2,mtemp3,vt)
941 .
' FLEXIBLE BODY INITIALIZATION FROM DMIG'/
942 .
' ---------------------- ')
9441100
FORMAT( /6x,
'FLEXIBLE BODY ID ',i10,1x,a
945 . /10x,
'NUMBER OF MODES ',i10
946 . /10x,
'NUMBER OF RIGID BODY MODES ',i10
947 . /10x,
'NUMBER OF NODES ',i10
948 . /10x,
'STABILITY TIME STEP ',1pe10.3)