Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

c_dabnew.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago
import madx-5.01.00
File size: 155.7 KB

Line
1	! The Full Polymorphic Package
2	! The module in this file is, to the best of our knowledge,
3	! the property of Lawrence Berkeley National Laboratory
4	! Its distribution and commercial usage may therefore be governed by the laws of the
5	! United States of America
6	module dabnew
7	use da_arrays
8	implicit none
9	! private
10	public
11	private daallno1,daall,damult,dasqrt,dacmut,dacma,DALINt,dafunt,dacctt
12	private dainvt,dapint,dadert,dacfuRt,dacfuIt,dacfut
13	private dadeb,dapac,dachk,damch,dadcd,dancd,hash,dehash
14	! public ppushstore,ppushprint,ppushGETN
15	! public DEALLOC_ALL,dainf
16	! PUBLIC DAINI,daall1
17	! PUBLIC DANOT,DAVAR,PPUSH,MTREE,DADAL,DACCT,DADER,DASQR
18	! PUBLIC DAMUL,DAADD,DACOP,DAINV,DAPIN,DAPEK,DAPOK,DACFU,DASUB,DACON,DADAL1
19	! PUBLIC DACLR,DACMU,DALIN,DAREA,DAPRI,DAABS,DAEPS,DATRA,MATINV,DAFUN
20	! PUBLIC DADIV,DADIC,DACDI,DACAD,DACSU,DASUC,DASHIFT,DARAN,DACFUR
21	! PUBLIC DACFUI,DAPRI77,DAREA77,GET_C_J,PPUSH1,DALLSTA,dacycle
22	! public count_da
23	! write_ removed
24	integer,private,parameter:: lsw=1
25	integer :: lda_max_used=0
26	! integer,private,parameter::nmax=400,lsw=1
27	! real(dp),private,parameter::tiny=c_1d_20
28	character(120),private :: line
29	contains
30	!******************************************************************************
31	! *
32	! *
33	! *
34	! DIFFERENTIAL ALGEBRA PACKAGE OF M. BERZ *
35	! **************************** *
36	! holy3 *
37	! *
38	! *
39	! *
40	! VERSION FOR MACHINE IN LINE THAT IS NOT COMMENTED OFF *
41	! TO CREATE DIFFERENT VERSIONS, USE THE PROGRAM 'VERSION' *
42	! *
43	! *
44	! *
45	! *
46	! THIS PACKAGE WAS INITIALLY WRITTEN BY PROF. M. BERZ WHILE AT *
47	! THE LAWRENCE BERKELEY LABORATORY. *
48	! IT HAS BEEN EXTENSIVELY MODIFIED BY THE MEMBERS OF THE ESG GROUP. *
49	! THEREFORE PROF. BERZ SHOULD NOT BE HELD RESPONSIBLE FOR ANY BUGS. *
50	! *
51	! NEW RULES OF THE GAME (EXHAUSTIVE) *
52	! ********************************** *
53	! THERE ARE NONE *
54	! *
55	!******************************************************************************
56	!
57	!
58	! THIS FILE CONTAINS ROUTINES TO PERFORM DIFFERENTIAL ALGEBRA (DA)
59	! AS AN OPTION, ALSO COMPONENTWISE ALGEBRA (CA) CAN BE PERFORMED.
60	! A DESCRIPTION OF THE INTERNAL ARRAYS USED BY THE ROUTINES CAN
61	! BE FOUND IN BLOCKDATA DABLD.
62	!
63	!
64	! SHORT REFERENCE CHART
65	! *********************
66	!
67	! THE PARAMETERS USED BELOW HAVE THE FOLLOWING MEANING:
68	!
69	! A,B: NAME OF INPUT DA VECTORS (INTEGER)
70	! C: NAME OF OUTPUT DA VECTOR (INTEGER)
71	! X,Y: NAME OF INPUT DA MATRIX (INTEGER(...))
72	! Z: NAME OF OUTPUT DA MATRIX (INTEGER(...))
73	!
74	! F: NAME OF A DA FUNCTION (CHARACTER(4))
75	! G: NAME OF EXTERNAL FUNCTION (real(dp))
76	! JJ: ARRAY OF EXPONENTS (INTEGER(20))
77	! O: ORDER (INTEGER)
78	! N: NUMBER OF VARIABLES (INTEGER)
79	! I,J,K: INTEGER NUMBER (INTEGER
80	! R,RA,RB: REAL NUMBERS (real(dp))
81	! H: ARRAY OF LENGTH LH (real(dp))
82	! U: OUTPUT UNIT NUMBER (INTEGER)
83	! T: COMMENT TEXT (CHARACTER(10))
84	!
85	!
86	! SUBROUTINES AND THEIR CALLING PARAMETERS
87	! ****************************************
88	!
89	! DAINI(O,N,U): INITIALIZES CONTROL ARRAYS AND SETS MAX. ORDER O AND
90	! MAX. NUMBER OF VARIABLES N. MUST BE CALLED BEFORE ANY
91	! OTHER DA ROUTINE CAN BE USED.
92	!
93	! DAALL(A,I,T,O,N): ALLOCATES SPACE FOR I VECTORS A. T: CHARACTER NAME
94	! DADAL(A,I): DEALLOCATES THE I VECTORS A.
95	!! DAVAR(A,R,I): MAKES A INDEPENDENT VARIABLE # I WITH INITIAL VALUE R
96	!! DACON(A,R): SETS A TO CONSTANT R
97	! DANOT(O): SETS NEW TRUNCATION ORDER O FOR DA OPERATIONS
98	! DAEPS(R): SETS NEW ZERO TOLERANCE EPSILON
99	!
100	!! DAPEK(A,JJ,R): RETURNS COEF R OF MONOMIAL WITH EXPONENTS JJ OF A
101	!! DAPOK(A,JJ,R): SETS COEF OF MONOMIAL WITH EXPONENTS JJ OF A TO R
102	!
103	!! DACOP(A,C): PERFORMS C = A
104	!! DAADD(A,B,C): PERFORMS C = A + B
105	!! DASUB(A,B,C): PERFORMS C = A - B
106	!! DAMUL(A,B,C): PERFORMS C = A * B
107	!! DADIV(A,B,C): PERFORMS C = A / B
108	!! DASQR(A,C): PERFORMS C = A^2 (SQUARE OF A)
109	!
110	!! DACAD(A,RA,C): PERFORMS C = A + RA
111	!! DACSU(A,RA,C): PERFORMS C = A - RA
112	!! DASUC(A,RA,C): PERFORMS C = RA - A
113	!! DACMU(A,RA,C): PERFORMS C = A * RA
114	!! DACDI(A,RA,C): PERFORMS C = A / RA
115	!! DADIC(A,RA,C): PERFORMS C = RA / A
116	!! DACMA(A,B,RB,C): PERFORMS C = A + RB*B
117	!!DAMULIN(A,B,RA,C,D,RB,C): PERFORMS C = ABRA + CDRB
118	!! DALIN(A,RA,B,RB,C): PERFORMS C = ARA + BRB
119	!! DAFUN(F,A,C): PERFORMS C = F(A) (DA FUNCTION)
120	!
121	!! DAABS(A,R): PERFORMS R = \|A\| (NORM OF A)
122	!! DACOM(A,B,R): PERFORMS R = \|A-B\| (NORM OF A-B)
123	!! DAPOS(A,C): PERFORMS C(I) = \|A(I)\| (MAKE SIGNS POSITIVE)
124	!
125	!! DACCT(X,I,Y,J,Z,K) CONCATENATES Z = X o Y; I,J,K: # OF VECTORS IN X,Y,
126	!! DAINV(X,I,Z,K) INVERTS Z = X^-1; I,J: # OF VECTORS IN X,Y
127	!! DAPIN(X,I,Z,K,JJ) PARTIALLY INVERTS Z = X^-1; I,J: # OF VECTORS IN X,Y,
128	! JJ: ARRAY; NONZERO ENTRIES DENOTE TO BE INVERTED LINES
129	!
130	!! DADER(I,A,C): PERFORMS C = DA/DI (DERIV. WITH RESPECT TO VARIABLE I
131	!! DAPOI(A,B,C,I): PERFORMS C = [A,B] (POISSON BRACKET, 2*I: # PHASEVARS
132	!! DACFU(A,G,C): MULTIPLIES COEFFICIENTS WITH FUNCTION G(JJ)
133	!
134	! DAPRI(A,U): PRINTS DA VECTOR A TO UNIT U
135	! DAREA(A,U): READS DA VECTOR A FROM UNIT U
136	! DADEB(U,T,I): DEBUGGER, DUMPS TO U. T: MEMO, I=0: RETURN, I=1:STOP
137	!! DARAN(A,R,seed): FILLS A WITH RANDOM NUMBERS. R: FILLFACTOR
138	! DANUM(O,N,I): COMPUTES NUMBER OF MONOMIALS IN N VAR THROUGH ORDER O
139	!
140	!
141	! ADDITIONAL ROUTINES THE USER DOES NOT NEED TO CALL:
142	!
143	! DAINF: RETURNS INFOS ABOUT A DA VECTOR PREVIOUSLY DECLARED
144	! DAPAC: PACKS DA VECTORS
145	! DACHK: CHECKS IF DA VECTORS HAVE COMPATIBLE ATTRIBUTES
146	! DCODE: TRANSFORMS DIGITS IN A CERTAIN BASE TO A DECIMAL INTEGER
147	! NCODE: EXTRACTS DIGITS IN A CERTAIN BASE FROM A DECIMAL INTEGER
148	!
149	!
150	! FURTHER WISHES
151	! **************
152	!
153	! - CHECK DAREA AND DAPRI FOR CA VECTORS
154	! - MAKE DAFUN USE DASQR
155	!
156	!
157	! BLOCKDATA DABLD
158	! ***************
159	!
160	!
161	! PARAMETERS:
162	!
163	! LDA: MAXIMUM NUMBER OF DA-VECTORS; CAN BE CHANGED QUITE ARBITRARILY
164	! LST: LENGTH OF MAIN STORAGE STACK; CAN BE CHANGED QUITE ARBITRARILY
165	! LEA: MAXIMUM NUMBER OF MONOMIALS; CAN BE INCREASED FOR LARGE NO,NV
166	! LIA: DIMENSION OF IA1,IA2; CAN BE INCREASED FOR LARGE NO,NV
167	! LNO: MAXIMUM ORDER; CAN BE INCREASED TO ABOUT 1000
168	! LNV: MAXIMUM NUMBER OF VARIABLES; CAN BE INCREASED TO ABOUT 1000
169	!
170	! ALL THE CHANGES IN THE VALUES OF PARAMETERS HAVE TO BE MADE BY GLOBAL
171	! SUBSTITUTIONS IN ALL SUBROUTINES.
172	!
173	! DANAME: NAME OF DA VECTOR
174	!
175	! CC: STACK OF DOUBLE PRECISON COEFFICIENTS
176	! i_1: FIRST CHARACTERISTIC INTEGER (CF DAINI)
177	! i_2: SECOND CHARACTERISTIC INTEGER (CF DAINI)
178	!
179	! IE1: CHARACTERISTIC INTEGER 1 OF UNPACKED REPRESENTATION (CF DAINI)
180	! IE2: CHARACTERISTIC INTEGER 2 OF UNPACKED REPRESENTATION (CF DAINI)
181	! IEO: ORDER OF ENTRY IN UNPACKED REPRESENTATION
182	! IA1: REVERSE TO IE1 (CF DAINI)
183	! IA2: REVERSE TO IE2 (CF DAINI)
184	!
185	! IDANO: ORDER OF DA VECTOR; IN CA, NUMBER OF COMPONENTS
186	! IDANV: NUMBER OF VARIABLES; IF 0, INDICATES CA VECTOR
187	! IDAPO: FIRST ADDRESS IN STACK
188	! IDALM: NUMBER OF RESERVED STACK POSITIONS
189	! IDALL: NUMBER OF MOMENTARILY REQUIRED STACK POSITIONS
190	!
191	! nda_dab: NUMBER OF DA VECTORS MOMENTARILY DEFINED
192	! NST: NUMBER OF STACK POSITIONS MOMENTARILY ALLOCATED
193	! NOMAX: MAXIMUM REQUESTED ORDER (CF DAINI)
194	! NVMAX: MAXIMUM REQUESTED NUMBER OF VARIABLES (CF DAINI)
195	! NMMAX: MAXIMUM NUMBER OF MONOMIALS FOR NOMAX, NVMAX (CF DAINI)
196	! NOCUT: MOMENTARY TRUNCATION ORDER
197	! EPS: TRUNCATION ACCURACY (CAN BE SET BY USER)
198	! EPSMAC: MANTISSA LENGTH OF MACHINE (PESSIMISTIC ESTIMATE)
199	!
200	!-----------------------------------------------------------------------------
201	!
202	subroutine change_package(i)
203	implicit none
204	integer, intent(in) :: i
205	if(i==2) then
206	lingyun_yang=.false.
207	elseif(i==1) then
208	lingyun_yang=.true.
209	else
210	write(6,*) " i = 1 or 2 "
211	write(6,*) " INPUT IGNORED "
212	endif
213	end subroutine change_package
214	subroutine daini(no,nv,iunit)
215	implicit none
216	! *****************************
217	!
218	! THIS SUBROUTINE SETS UP THE MAJOR ORDERING AND ADDRESSING ARRAYS IN
219	! COMMON BLOCK DAINI. IF IUNIT > 0, THE ARRAYS WILL BE PRINTED TO UNIT
220	! NUMBER IUNIT. AN EXAMPLE FOR THE ARRAYS GENERATED BY DAINI CAN BE
221	! FOUND AFTER THE ROUTINE.
222	!
223	!-----------------------------------------------------------------------------
224	!-----------------------------------------------------------------------------
225	! COMMON / DASCR / IS(20), RS(20)
226	!integer idao,is,iscrri
227	!real(dp) rs
228	!common/dascr/is(100),rs(100),iscrri(100),idao
229	!-----------------------------------------------------------------------------
230	!
231	integer i,iall,ibase,ic1,ic2,icmax,io1,io2,iout,iunit,jd,jjj,jjjj,jl,js,&
232	nn,no,nv,ipause,mypauses
233	integer,dimension(lnv+1)::n
234	integer,dimension(0:lnv)::k
235	integer,dimension(lnv)::j,jj
236	character(10) aa
237	!
238	!frs if(eps.le.zero) eps=c_1d_38
239	! if(EPS.le.zero) eps=c_1d_90
240	!frs epsmac=c_1d_7
241	if((.not.C_%STABLE_DA)) then
242	if(c_%watch_user) then
243	write(6,*) "big problem in dabnew ", sqrt(crash)
244	endif
245	return
246	endif
247	last_tpsa=2
248	if(nv.eq.0) return
249	call alloc_all(no,nv)
250	ndamaxi=0
251	!
252	do i=1, lda
253	allvec(i) = .false.
254	enddo
255	nhole=0
256	!*****************************************
257	! INITIALIZING VARIABLES IN COMMON / DA /
258	! ***************************************
259	!
260	nda_dab = 0
261	nst0 = 0
262	nomax = no
263	nvmax = nv
264	call danum(no,nv,nmmax)
265	nocut = no
266	lfi = 0
267	!
268	do i=0,lia
269	ia1(i) = 0
270	ia2(i) = 0
271	enddo
272	!
273	! do i=1,100
274	! is(i) = 0
275	! enddo
276	!
277	if(nv.gt.lnv.or.no.gt.lno) then
278	write(6,*) 'ERROR IN SUBROUTINE DAINI, NO, NV = ',no,nv
279	ipause=mypauses(1,line)
280	call dadeb !(31,'ERR DAINI ',1)
281	endif
282	!
283	ibase = no+1
284	js = nv/2
285	if(float(ibase)**((nv+1)/2).gt.float(lia)) then
286	write(line,'(a12,i4,a7,i4,a21,i4)') 'ERROR, NO = ',no,', NV = ',nv,' TOO LARGE FOR LIA = ',lia
287	ipause=mypauses(2,line)
288	call dadeb !(31,'ERR DAINI ',1)
289	endif
290	!
291	icmax = 0
292	nn = 0
293	k(0) = 0
294	!
295	do io2=0,no
296	! ***************
297	!
298	n(1) = io2
299	jl = 0
300	jd = 1
301	!
302	50 jl = jl + jd
303	!
304	!old
305	! IF(JL.EQ.0) THEN
306	!old
307	! modified according to Wu Ying
308	!
309	if(jl.le.0) then
310	!
311	goto 100
312	elseif(jd.eq.1) then
313	j(jl) = 0
314	else
315	j(jl) = j(jl) + 1
316	endif
317	!
318	k(jl) = k(jl-1)*ibase + j(jl)
319	n(jl+1) = n(jl) - j(jl)
320	!
321	if(j(jl).gt.n(jl)) then
322	jd = -1
323	goto 50
324	elseif(jl.lt.js) then
325	jd = 1
326	goto 50
327	else
328	j(jl) = n(jl)
329	k(jl) = k(jl-1)*ibase + j(jl)
330	ic2 = k(jl)
331	icmax = max(icmax,ic2)
332	k(jl) = 0
333	!
334	ia2(ic2) = nn
335	!
336	do io1=0,no-io2
337	! ******************
338	!
339	n(js+1) = io1
340	jd = 1
341	!
342	70 jl = jl + jd
343	!
344	if(jl.eq.js) then
345	goto 80
346	elseif(jd.eq.1) then
347	j(jl) = 0
348	else
349	j(jl) = j(jl) + 1
350	endif
351	!
352	k(jl) = k(jl-1)*ibase + j(jl)
353	n(jl+1) = n(jl) - j(jl)
354	!
355	if(j(jl).gt.n(jl)) then
356	jd = -1
357	goto 70
358	elseif(jl.lt.nv) then
359	jd = 1
360	goto 70
361	else
362	jd = -1
363	j(jl) = n(jl)
364	k(jl) = k(jl-1)*ibase + j(jl)
365	ic1 = k(jl)
366	icmax = max(icmax,ic1)
367	nn = nn + 1
368	!
369	if(etiennefix) then
370	if(nn<=lea) then ! Etienne
371	ie1(nn) = ic1
372	ie2(nn) = ic2
373	endif ! Etienne
374	if(nn<=lst) then ! Etienne
375	i_1 (nn) = ic1
376	i_2 (nn) = ic2
377	endif
378	if(ic2.eq.0) ia1(ic1) = nn
379	if(nn<=lea) then ! Etienne
380	ieo(nn) = io1 + io2
381	endif ! Etienne
382	!
383	else
384	ie1(nn) = ic1
385	ie2(nn) = ic2
386	i_1 (nn) = ic1
387	i_2 (nn) = ic2
388	if(ic2.eq.0) ia1(ic1) = nn
389	ieo(nn) = io1 + io2
390	!
391	endif
392	goto 70
393	endif
394	!
395	80 continue
396	enddo
397	!
398	jd = -1
399	goto 50
400	endif
401	!
402	100 continue
403	enddo
404	!
405	if(nn.gt.lea.and.(.not.etiennefix)) then
406	write(line,'(a21,i4,a12)') 'ERROR IN DAINI, NN = ',nn,' EXCEEDS LEA'
407	ipause=mypauses(3,line)
408	call dadeb !(31,'ERR DAINI ',1)
409	endif
410	!
411	! ALLOCATING SCRATCH VARIABLES
412	! ****************************
413	!
414	iall = 0
415	call daall1(iall,'$$UNPACK$$',nomax,nvmax)
416	!
417	do i=0,nomax
418	aa = '$$MUL $$'
419	write(aa(6:10),'(I5)') i
420	iall = 0
421	! CALL DAALL(IALL,1,AA,I,NVMAX)
422	call daall1(iall,aa,nomax,nvmax)
423	enddo
424	!
425	idall(1) = nmmax
426	!
427	! DOUBLE CHECKING ARRAYS IE1,IE2,IA1,IA2
428	! **************************************
429	!
430	do i=1,nmmax
431	!
432	jjj = ia1(ie1(i)) + ia2(ie2(i))
433	if(jjj.ne.i) then
434	write(line,'(a48,i4)') 'ERROR IN DAINI IN ARRAYS IE1,IE2,IA1,IA2 AT I = ',i
435	ipause=mypauses(4,line)
436	call dadeb !(31,'ERR DAINI ',1)
437	endif
438	!
439	enddo
440	!
441	if(iunit.eq.0) return
442	!
443	write(line,'(a32)') 'ARRAY SETUP DONE, BEGIN PRINTING'
444	ipause=mypauses(5,line)
445	!
446	iout = 32
447	open(iout,file='DAINI.DAT',status='NEW')
448	!CRAY OPEN(IOUT,FILE='DAINI',STATUS='UNKNOWN',FORM='FORMATTED') *CRAY
449	!CRAY REWIND IOUT *CRAY
450	!
451	write(iout,'(/A/A/)') ' ARRAYS i_1 THROUGH i_20, IE1,IE2,IEO **********************************'
452	do i=1,nmmax
453	call dancd(ie1(i),ie2(i),jj)
454	write(iout,'(1X,I5,2X,4(5i2,1X),3I6)') i,(jj(jjjj),jjjj=1,lnv),ie1(i),ie2(i),ieo(i)
455	enddo
456	!
457	write(iout,'(/A/A/)') ' ARRAYS IA1,IA2 **************'
458	do i=0,icmax
459	write(iout,'(3i10)') i,ia1(i),ia2(i)
460	enddo
461	!
462	return
463	end subroutine daini
464	! subroutine daexter
465	! implicit none
466	! *****************************
467	!
468	!-----------------------------------------------------------------------------
469	! integer i_1,i_2,ia1,ia2,idall,idalm,idano,idanv,idapo,ie1,ie2,ieo,ifi,lfi,nda,ndamaxi,nmmax,nocut,nomax,nst,nvmax
470	!-----------------------------------------------------------------------------
471	!
472	! integer i
473	!
474	! if(.not.notallocated) then
475	! do i=1, lda
476	! allvec(i)=.false.
477	! enddo
478	! endif
479	! return
480	! end subroutine daexter
481	subroutine dallsta(ldanow)
482	implicit none
483	! *****************************
484	!
485	!-----------------------------------------------------------------------------
486	!
487	integer i,ldanow
488	!
489	ldanow=0
490	do i=1, lda
491	if(allvec(i)) ldanow=ldanow+1
492	enddo
493	w_p=0
494	w_p%nc=1
495	write(w_p%c(1),'(a11,i8)') ' ALLOCATED ',ldanow
496	w_p%fc='(1((1X,A72)))'
497	!CALL WRITE_i
498	return
499	end subroutine dallsta
500	!
501	! EXAMPLE: ARRAYS i_1 THROUGH i_20, IE1,IE2,IEO (NOMAX=3,NVMAX=4)
502	! *************************************************************
503	! I i_1 THROUGH i_20 IE1 IE2 IEO
504	! 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
505	! 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
506	! 3 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 1
507	! 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 2
508	! 5 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 2
509	! 6 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 2
510	! 7 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 3
511	! 8 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 3
512	! 9 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9 0 3
513	! 10 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 3
514	! 11 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1
515	! 12 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 2
516	! 13 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 4 1 2
517	! 14 2 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 2 1 3
518	! 15 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 5 1 3
519	! 16 0 2 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 8 1 3
520	! 17 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 4 1
521	! 18 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 4 2
522	! 19 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 4 4 2
523	! 20 2 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 2 4 3
524	! 21 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 5 4 3
525	! 22 0 2 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 8 4 3
526	! 23 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 2 2
527	! 24 1 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 1 2 3
528	! 25 0 1 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 4 2 3
529	! 26 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 5 2
530	! 27 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 5 3
531	! 28 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 4 5 3
532	! 29 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 8 2
533	! 30 1 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 1 8 3
534	! 31 0 1 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 4 8 3
535	! 32 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 3 3
536	! 33 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 0 6 3
537	! 34 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 0 0 9 3
538	! 35 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 12 3
539	!
540	! ARRAYS IA1,IA2
541	! **************
542	! I IA1 IA2
543	! 0 1 0 IE1,IE2 AND IA1,IA2 ALLOW THE EASY COMPUTATION
544	! 1 2 10 OF THE ADDRESS OF THE PRODUCT OF TWO MONOMIALS.
545	! 2 4 22 LET IX AND IY BE THE POSITIONS OF THE TWO
546	! 3 7 31 FACTORS. THEN THE POSITION IZ OF THE PRODUCT OF
547	! 4 3 16 THE TWO FACTORS IS GIVEN BY
548	! 5 5 25
549	! 6 8 32 IZ = IA1(IE1(IX)+IE1(IY)) + IA2(IE2(IX)+IE2(IY))
550	! 7 0 0
551	! 8 6 28
552	! 9 9 33 THE OTHER VARIABLES SET BY DAINI WOULD HAVE THE
553	! 10 0 0 VALUES
554	! 11 0 0
555	! 12 10 34 NOMAX = 3, NVMAX = 4, NMMAX = 35
556	!
557	subroutine daallno1(ic,ccc)
558	implicit none
559	! ********************************
560	!
561	! THIS SUBROUTINE ALLOCATES STORAGE FOR A DA VECTOR WITH
562	! ORDER NOmax AND NUMBER OF VARIABLES NVmax
563	!
564	!-----------------------------------------------------------------------------
565	!
566	logical(lp) incnda
567	integer ind,ndanum,no,nv,ic,ipause,mypauses
568	character(10) c,ccc
569	! if((.not.C_%STABLE_DA)) then
570	! if(c_%watch_user) then
571	! write(6,*) "big problem in dabnew ", sqrt(crash)
572	! endif
573	! return
574	! endif
575	!
576	no=nomax
577	nv=nvmax
578	ind = 1
579	if(ic.gt.0.and.ic.le.nda_dab) then
580	! DANAME(IC(I)) = C
581	! IF(IDANO(IC(I)).EQ.NO.AND.IDANV(IC(I)).EQ.NV) THEN
582	else
583	if(nv.ne.0.and.(no.gt.nomax.or.nv.gt.nvmax)) then
584	write(line,'(a23,i4,a14,i4,1x,i4,a16,i4,1x,i4)') 'ERROR IN DAALL, VECTOR ',c,' HAS NO, NV = ',no,nv, &
585	&' NOMAX, NVMAX = ',nomax,nvmax
586	ipause=mypauses(7,line)
587	call dadeb !(31,'ERR DAALL ',1)
588	endif
589	!
590	if(nhole.gt.0) then
591	ind=nda_dab
592	20 if (allvec(ind)) then
593	ind = ind - 1
594	goto 20
595	endif
596	incnda = .false.
597	nhole=nhole-1
598	else
599	incnda = .true.
600	nda_dab = nda_dab + 1
601	ind=nda_dab
602	if(nda_dab.gt.lda) then
603	write(line,'(a50)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED'
604	ipause=mypauses(8,line)
605	call dadeb !(31,'ERR DAALL ',1)
606	endif
607	endif
608	if(ind>lda_max_used) lda_max_used=ind
609	if(ind>lda) then
610	write(6,*) "ind>lda ",lda,ind
611	print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA
612	stop
613	endif
614	allvec(ind) = .true.
615	ic = ind
616	!
617	if(nv.ne.0) then
618	call danum(no,nv,ndanum)
619	else
620	ndanum = no
621	endif
622	c = ccc
623	write(c(6:10),'(I5)') 1
624	daname(ind) = c
625	if (incnda) then
626	if(ind.gt.nomax+2) then
627	idano(ind) = nomax
628	idanv(ind) = nvmax
629	idapo(ind) = nst0 + 1
630	idalm(ind) = nmmax
631	idall(ind) = 0
632	nst0 = nst0 + nmmax
633	else
634	idano(ind) = no
635	idanv(ind) = nv
636	idapo(ind) = nst0 + 1
637	idalm(ind) = ndanum
638	idall(ind) = 0
639	nst0 = nst0 + ndanum
640	endif
641	endif
642	!
643	if(nst0.gt.lst) then
644	w_p=0
645	w_p%nc=5
646	w_p%fc='(4(1X,a72,/),(1X,a72))'
647	w_p%c(1)= 'ERROR IN DAALL, STACK EXHAUSTED '
648	w_p%c(2)= ' NST0,LST '
649	write(w_p%c(3),'(i8,1x,i8)') NST0,LST
650	w_p%c(4)= ' NDA,NDANUM,NDA*NDANUM '
651	write(w_p%c(5),'(i8,1x,i8,1x,i8)') nda_dab,ndanum,nda_dab*ndanum
652	!CALL WRITE_E(125)(125)
653	call dadeb !(31,'ERR DAALL ',1)
654	endif
655	!
656	if(nv.eq.0.or.nomax.eq.1) then
657	call daclr(ic)
658	idall(ic) = idalm(ic)
659	endif
660	endif
661	!
662	if(nda_dab.gt.ndamaxi) ndamaxi=nda_dab
663	return
664	end subroutine daallno1
665	subroutine daall(ic,l,ccc,no,nv)
666	implicit none
667	! ********************************
668	!
669	! THIS SUBROUTINE ALLOCATES STORAGE FOR A DA VECTOR WITH
670	! ORDER NO AND NUMBER OF VARIABLES NV
671	!
672	!-----------------------------------------------------------------------------
673	!
674	logical(lp) incnda
675	integer i,ind,l,ndanum,no,nv,ipause,mypauses
676	integer,dimension(:)::ic
677	character(10) c,ccc
678	! if((.not.C_%STABLE_DA)) then
679	! if(c_%watch_user) then
680	! write(6,*) "big problem in dabnew ", sqrt(crash)
681	! endif
682	! return
683	! endif
684	!
685	ind = 1
686	do i=1,l
687	if(ic(i).gt.0.and.ic(i).le.nda_dab) then
688	! DANAME(IC(I)) = C
689	! IF(IDANO(IC(I)).EQ.NO.AND.IDANV(IC(I)).EQ.NV) THEN
690	else
691	if(nv.ne.0.and.(no.gt.nomax.or.nv.gt.nvmax)) then
692	write(line,'(a23,i4,a14,i4,1x,i4,a16,i4,1x,i4)') 'ERROR IN DAALL, VECTOR ',c,' HAS NO, NV = ',no,nv, &
693	&' NOMAX, NVMAX = ',nomax,nvmax
694	ipause=mypauses(9,line)
695	call dadeb !(31,'ERR DAALL ',1)
696	endif
697	!
698	if(nhole.gt.0) then
699	ind=nda_dab
700	20 if (allvec(ind)) then
701	ind = ind - 1
702	goto 20
703	endif
704	incnda = .false.
705	nhole=nhole-1
706	else
707	incnda = .true.
708	nda_dab = nda_dab + 1
709	ind=nda_dab
710	if(nda_dab.gt.lda) then
711	write(6,'(a50)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED'
712	! ipause=mypauses(10,line)
713	call dadeb !(31,'ERR DAALL ',1)
714	stop 111
715	endif
716	endif
717	!write(30,*) no,ind,lda,size(allvec)
718	if(ind>lda_max_used) lda_max_used=ind
719	if(ind>lda) then
720	write(6,*) "ind>lda ",lda,ind
721	print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA
722	stop
723	endif
724	allvec(ind) = .true.
725	ic(i) = ind
726	!
727	if(nv.ne.0) then
728	call danum(no,nv,ndanum)
729	else
730	ndanum = no
731	endif
732	c = ccc
733	if(l.ne.1) write(c(6:10),'(I5)') i
734	daname(ind) = c
735	if (incnda) then
736	if(ind.gt.nomax+2) then
737	idano(ind) = nomax
738	idanv(ind) = nvmax
739	idapo(ind) = nst0 + 1
740	idalm(ind) = nmmax
741	idall(ind) = 0
742	nst0 = nst0 + nmmax
743	else
744	idano(ind) = no
745	idanv(ind) = nv
746	idapo(ind) = nst0 + 1
747	idalm(ind) = ndanum
748	idall(ind) = 0
749	nst0 = nst0 + ndanum
750	endif
751	endif
752	!
753	if(nst0.gt.lst) then
754	w_p=0
755	w_p%nc=5
756	w_p%fc='(4(1X,a72,/),(1X,a72))'
757	w_p%c(1)= 'ERROR IN DAALL, STACK EXHAUSTED '
758	w_p%c(2)= ' NST,LST '
759	write(w_p%c(3),'(i8,1x,i8)') NST0,LST
760	w_p%c(4)= ' NDA,NDANUM,NDA*NDANUM '
761	write(w_p%c(5),'(i8,1x,i8,1x,i8)') nda_dab,ndanum,nda_dab*ndanum
762	!CALL WRITE_E(125)(126)
763	call dadeb !(31,'ERR DAALL ',1)
764	endif
765	!
766	! IF(NV.EQ.0) THEN
767	if(nv.eq.0.or.nomax.eq.1) then
768	call daclr(ic(i))
769	idall(ic(i)) = idalm(ic(i))
770	endif
771	endif
772	enddo
773	!
774	if(nda_dab.gt.ndamaxi) ndamaxi=nda_dab
775	return
776	end subroutine daall
777	subroutine daall1(ic,ccc,no,nv)
778	implicit none
779	! ********************************
780	!
781	! THIS SUBROUTINE ALLOCATES STORAGE FOR A DA VECTOR WITH
782	! ORDER NO AND NUMBER OF VARIABLES NV
783	!
784	!-----------------------------------------------------------------------------
785	!
786	logical(lp) incnda
787	integer ic,ind,ndanum,no,nv,ipause,mypauses
788	character(10) c,ccc
789	! if((.not.C_%STABLE_DA)) then
790	! if(c_%watch_user) then
791	! write(6,*) "big problem in dabnew ", sqrt(crash)
792	! endif
793	! return
794	! endif
795	!
796	ind = 1
797	if(ic.gt.0.and.ic.le.nda_dab) then
798	! DANAME(ic) = C
799	! IF(IDANO(ic).EQ.NO.AND.IDANV(ic).EQ.NV) THEN
800	else
801	if(nv.ne.0.and.(no.gt.nomax.or.nv.gt.nvmax)) then
802	write(line,'(a23,i4,a14,i4,1x,i4,a16,i4,1x,i4)') 'ERROR IN DAALL, VECTOR ',c,' HAS NO, NV = ',no,nv, &
803	&' NOMAX, NVMAX = ',nomax,nvmax
804	ipause=mypauses(11,line)
805	call dadeb !(31,'ERR DAALL ',1)
806	endif
807	!
808	if(nhole.gt.0) then
809	ind=nda_dab
810	20 if (allvec(ind)) then
811	ind = ind - 1
812	goto 20
813	endif
814	incnda = .false.
815	nhole=nhole-1
816	else
817	incnda = .true.
818	nda_dab = nda_dab + 1
819	ind=nda_dab
820	if(nda_dab.gt.lda) then
821	write(line,'(a50)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED'
822	ipause=mypauses(12,line)
823	call dadeb !(31,'ERR DAALL ',1)
824	endif
825	endif
826	if(ind>lda_max_used) lda_max_used=ind
827	if(ind>lda) then
828	write(6,*) "ind>lda ",lda,ind
829	print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA
830	stop
831	endif
832	allvec(ind) = .true.
833	ic = ind
834	!
835	if(nv.ne.0) then
836	call danum(no,nv,ndanum)
837	else
838	ndanum = no
839	endif
840	c = ccc
841	write(c(6:10),'(I5)') 1
842	daname(ind) = c
843	if (incnda) then
844	if(ind.gt.nomax+2) then
845	idano(ind) = nomax
846	idanv(ind) = nvmax
847	idapo(ind) = nst0 + 1
848	idalm(ind) = nmmax
849	idall(ind) = 0
850	nst0 = nst0 + nmmax
851	else
852	idano(ind) = no
853	idanv(ind) = nv
854	idapo(ind) = nst0 + 1
855	idalm(ind) = ndanum
856	idall(ind) = 0
857	nst0 = nst0 + ndanum
858	endif
859	endif
860	!
861	if(nst0.gt.lst) then
862	w_p=0
863	w_p%nc=5
864	w_p%fc='(4(1X,a72,/),(1X,a72))'
865	w_p%c(1)= 'ERROR IN DAALL, STACK EXHAUSTED '
866	w_p%c(2)= ' NST,LST '
867	write(w_p%c(3),'(i8,1x,i8)') NST0,LST
868	w_p%c(4)= ' NDA,NDANUM,NDA*NDANUM '
869	write(w_p%c(5),'(i8,1x,i8,1x,i8)') nda_dab,ndanum,nda_dab*ndanum
870	!CALL WRITE_E(125)(127)
871	call dadeb !(31,'ERR DAALL ',1)
872	endif
873	!
874	! IF(NV.EQ.0) THEN
875	if(nv.eq.0.or.nomax.eq.1) then
876	call daclr(ic)
877	idall(ic) = idalm(ic)
878	endif
879	endif
880	!
881	if(nda_dab.gt.ndamaxi) ndamaxi=nda_dab
882	return
883	end subroutine daall1
884	subroutine daall0(ic)
885	implicit none
886	! ********************************
887	!
888	! THIS SUBROUTINE ALLOCATES STORAGE FOR A DA VECTOR WITH
889	! ORDER NO AND NUMBER OF VARIABLES NV
890	!
891	!-----------------------------------------------------------------------------
892	!
893	logical(lp) incnda
894	integer ic,ind,ndanum,no,nv,ipause,mypauses
895	character(10) c,ccc
896	ccc=' '
897	no=nomax
898	nv=nvmax
899	! if((.not.C_%STABLE_DA)) then
900	! if(c_%watch_user) then
901	! write(6,*) "big problem in dabnew ", sqrt(crash)
902	! endif
903	! return
904	! endif
905	!
906	ind = 1
907	if(ic.gt.0.and.ic.le.nda_dab) then
908	! DANAME(ic) = C
909	! IF(IDANO(ic).EQ.NO.AND.IDANV(ic).EQ.NV) THEN
910	else
911	if(nv.ne.0.and.(no.gt.nomax.or.nv.gt.nvmax)) then
912	write(line,'(a23,i4,a14,i4,1x,i4,a16,i4,1x,i4)') 'ERROR IN DAALL, VECTOR ',c,' HAS NO, NV = ',no,nv, &
913	&' NOMAX, NVMAX = ',nomax,nvmax
914	ipause=mypauses(11,line)
915	call dadeb !(31,'ERR DAALL ',1)
916	endif
917	!
918	if(nhole.gt.0) then
919	ind=nda_dab
920	20 if (allvec(ind)) then
921	ind = ind - 1
922	goto 20
923	endif
924	incnda = .false.
925	nhole=nhole-1
926	else
927	incnda = .true.
928	nda_dab = nda_dab + 1
929	ind=nda_dab
930	if(nda_dab.gt.lda) then
931	write(line,'(a50)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED'
932	ipause=mypauses(12,line)
933	call dadeb !(31,'ERR DAALL ',1)
934	endif
935	endif
936	if(ind>lda_max_used) lda_max_used=ind
937	if(ind>lda) then
938	write(6,*) "ind>lda ",lda,ind
939	print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA
940	stop
941	endif
942	allvec(ind) = .true.
943	ic = ind
944	!
945	if(nv.ne.0) then
946	call danum(no,nv,ndanum)
947	else
948	ndanum = no
949	endif
950	c = ccc
951	write(c(6:10),'(I5)') 1
952	daname(ind) = c
953	if (incnda) then
954	if(ind.gt.nomax+2) then
955	idano(ind) = nomax
956	idanv(ind) = nvmax
957	idapo(ind) = nst0 + 1
958	idalm(ind) = nmmax
959	idall(ind) = 0
960	nst0 = nst0 + nmmax
961	else
962	idano(ind) = no
963	idanv(ind) = nv
964	idapo(ind) = nst0 + 1
965	idalm(ind) = ndanum
966	idall(ind) = 0
967	nst0 = nst0 + ndanum
968	endif
969	endif
970	!
971	if(nst0.gt.lst) then
972	w_p=0
973	w_p%nc=5
974	w_p%fc='(4(1X,a72,/),(1X,a72))'
975	w_p%c(1)= 'ERROR IN DAALL, STACK EXHAUSTED '
976	w_p%c(2)= ' NST,LST '
977	write(w_p%c(3),'(i8,1x,i8)') NST0,LST
978	w_p%c(4)= ' NDA,NDANUM,NDA*NDANUM '
979	write(w_p%c(5),'(i8,1x,i8,1x,i8)') nda_dab,ndanum,nda_dab*ndanum
980	!CALL WRITE_E(125)(127)
981	call dadeb !(31,'ERR DAALL ',1)
982	endif
983	!
984	! IF(NV.EQ.0) THEN
985	if(nv.eq.0.or.nomax.eq.1) then
986	call daclr(ic)
987	idall(ic) = idalm(ic)
988	endif
989	endif
990	!
991	if(nda_dab.gt.ndamaxi) ndamaxi=nda_dab
992	return
993	end subroutine daall0
994	!
995	subroutine dadal(idal,l)
996	implicit none
997	! ************************
998	!
999	! THIS SUBROUTINE DEALLOCATES THE VECTORS IDAL
1000	!
1001	!-----------------------------------------------------------------------------
1002	!
1003	integer i,l,ipause,mypauses
1004	integer,dimension(:)::idal
1005	!
1006	! if((.not.C_%STABLE_DA)) then
1007	! if(c_%watch_user) then
1008	! write(6,*) "big problem in dabnew ", sqrt(crash)
1009	! endif
1010	! return
1011	! endif
1012	do i=l,1,-1
1013	if(idal(i).le.nomax+2.or.idal(i).gt.nda_dab) then
1014	write(line,'(a38,i8,1x,i8)') 'ERROR IN ROUTINE DADAL, IDAL(I),NDA = ',idal(i),nda_dab
1015	ipause=mypauses(13,line)
1016	call dadeb !(31,'ERR DADAL ',1)
1017	endif
1018	if(idal(i).eq.nda_dab) then
1019	! deallocate
1020	nst0 = idapo(nda_dab) - 1
1021	nda_dab = nda_dab - 1
1022	else
1023	nhole=nhole+1
1024	endif
1025	allvec(idal(i)) = .false.
1026	! IDANO(IDAL(I)) = 0
1027	! IDANV(IDAL(I)) = 0
1028	! IDAPO(IDAL(I)) = 0
1029	! IDALM(IDAL(I)) = 0
1030	idall(idal(i)) = 0
1031	idal(i) = 0
1032	enddo
1033	return
1034	end subroutine dadal
1035	subroutine dadal1(idal)
1036	implicit none
1037	! ************************
1038	!
1039	! THIS SUBROUTINE DEALLOCATES THE VECTORS IDAL
1040	!
1041	!-----------------------------------------------------------------------------
1042	!
1043	integer idal,ipause,mypauses
1044	!
1045	! if((.not.C_%STABLE_DA)) then
1046	! if(c_%watch_user) then
1047	! write(6,*) "big problem in dabnew ", sqrt(crash)
1048	! endif
1049	! return
1050	! endif
1051	if(idal.le.nomax+2.or.idal.gt.nda_dab) then
1052	write(line,'(a35,i8,1x,i8)') 'ERROR IN ROUTINE DADAL, IDAL,NDA = ',idal,nda_dab
1053	ipause=mypauses(14,line)
1054	call dadeb !(31,'ERR DADAL ',1)
1055	endif
1056	if(idal.eq.nda_dab) then
1057	! deallocate
1058	nst0 = idapo(nda_dab) - 1
1059	nda_dab = nda_dab - 1
1060	else
1061	nhole=nhole+1
1062	endif
1063	allvec(idal) = .false.
1064	! IDANO(IDAL(I)) = 0
1065	! IDANV(IDAL(I)) = 0
1066	! IDAPO(IDAL(I)) = 0
1067	! IDALM(IDAL(I)) = 0
1068	idall(idal) = 0
1069	idal = 0
1070	return
1071	end subroutine dadal1
1072	subroutine count_da(n)
1073	implicit none
1074	! ************************
1075	!
1076	! THIS SUBROUTINE counts allocate da
1077	!
1078	!-----------------------------------------------------------------------------
1079	!
1080	integer i,n
1081	!
1082	n=0
1083	do i=1,lda
1084	if(allvec(i)) n=n+1
1085	enddo
1086	return
1087	end subroutine count_da
1088	subroutine davar(ina,ckon,i)
1089	implicit none
1090	! ****************************
1091	!
1092	! THIS SUBROUTINE DECLARES THE DA VECTOR
1093	! AS THE INDEPENDENT VARIABLE NUMBER I.
1094	!
1095	!-----------------------------------------------------------------------------
1096	!
1097	integer i,ibase,ic1,ic2,illa,ilma,ina,inoa,inva,ipoa,ipause,mypauses
1098	real(dp) ckon
1099	!
1100	if((.not.C_%STABLE_DA)) then
1101	if(c_%watch_user) then
1102	write(6,*) "big problem in dabnew ", sqrt(crash)
1103	endif
1104	return
1105	endif
1106	call dainf(ina,inoa,inva,ipoa,ilma,illa)
1107	if((.not.C_%STABLE_DA)) then
1108	if(c_%watch_user) then
1109	write(6,*) "big problem in dabnew ", sqrt(crash)
1110	endif
1111	return
1112	endif
1113	!
1114	!
1115	if(i.gt.inva) then
1116	write(line,'(a20,i8,a16,i8)') 'ERROR IN DAVAR, I = ',i,' EXCEEDS INVA = ',inva
1117	ipause=mypauses(14,line)
1118	call dadeb !(31,'ERR DAVAR ',1)
1119	endif
1120	!
1121	if(nomax.eq.1) then
1122	if(i.gt.inva) then
1123	print*,'ERROR IN DAVAR, I = ',i,' EXCEEDS INVA = ',inva
1124	! call dadeb !(31,'ERR DAVAR3',1)
1125	endif
1126	call daclr(ina)
1127	cc(ipoa) = ckon
1128	cc(ipoa+i) = 1.0_dp
1129	return
1130	endif
1131	ibase = nomax+1
1132	!
1133	if(i.gt.(nvmax+1)/2) then
1134	ic1 = 0
1135	ic2 = ibase**(i-(nvmax+1)/2-1)
1136	else
1137	ic1 = ibase**(i-1)
1138	ic2 = 0
1139	endif
1140	!
1141	if(abs(ckon).gt.eps) then
1142	idall(ina) = 2
1143	cc(ipoa) = ckon
1144	i_1(ipoa) = 0
1145	i_2(ipoa) = 0
1146	!
1147	cc(ipoa+1) = 1.0_dp
1148	i_1(ipoa+1) = ic1
1149	i_2(ipoa+1) = ic2
1150	else
1151	idall(ina) = 1
1152	cc(ipoa) = 1.0_dp
1153	i_1(ipoa) = ic1
1154	i_2(ipoa) = ic2
1155	endif
1156	!
1157	return
1158	end subroutine davar
1159	!
1160	subroutine dacon(ina,ckon)
1161	implicit none
1162	! **************************
1163	!
1164	! THIS SUBROUTINE SETS THE VECTOR C TO THE CONSTANT CKON
1165	!
1166	!-----------------------------------------------------------------------------
1167	!
1168	integer illa,ilma,ina,inoa,inva,ipoa
1169	real(dp) ckon
1170	!
1171	if((.not.C_%STABLE_DA)) then
1172	if(c_%watch_user) then
1173	write(6,*) "big problem in dabnew ", sqrt(crash)
1174	endif
1175	return
1176	endif
1177	call dainf(ina,inoa,inva,ipoa,ilma,illa)
1178	if((.not.C_%STABLE_DA)) then
1179	if(c_%watch_user) then
1180	write(6,*) "big problem in dabnew ", sqrt(crash)
1181	endif
1182	return
1183	endif
1184	!
1185	!
1186	if(nomax.eq.1) then
1187	call daclr(ina)
1188	cc(ipoa) = ckon
1189	return
1190	endif
1191	idall(ina) = 1
1192	cc(ipoa) = ckon
1193	i_1(ipoa) = 0
1194	i_2(ipoa) = 0
1195	if(abs(ckon).lt.eps) idall(ina) = 0
1196	!
1197	return
1198	end subroutine dacon
1199	!
1200	subroutine danot(not)
1201	implicit none
1202	! *********************
1203	!
1204	! THIS SUBROUTINE RESETS THE TRUNCATION ORDER NOCUT TO A NEW VALUE
1205	!
1206	!-----------------------------------------------------------------------------
1207	!
1208	integer not,ipause,mypauses
1209	!
1210	if((.not.C_%STABLE_DA)) then
1211	if(c_%watch_user) then
1212	write(6,*) "big problem in dabnew ", sqrt(crash)
1213	endif
1214	return
1215	endif
1216	if(not.gt.nomax) then
1217	write(line,'(a15,i8,a17,i8)') 'ERROR, NOCUT = ',nocut,' EXCEEDS NOMAX = ',nomax
1218	ipause=mypauses(15,line)
1219	call dadeb !(31,'ERR DANOT ',1)
1220	endif
1221	!
1222	nocut = not
1223	!
1224	return
1225	end subroutine danot
1226	! subroutine getdanot(not)
1227	! implicit none
1228	! *********************
1229	!
1230	! THIS SUBROUTINE RESETS THE TRUNCATION ORDER NOCUT TO A NEW VALUE
1231	!
1232	!-----------------------------------------------------------------------------
1233	!
1234	! integer not,ipause,mypauses
1235	!
1236	! if(not.gt.nomax) then
1237	! write(line,'(a15,i8,a17,i8)') 'ERROR, NOCUT = ',nocut,' EXCEEDS NOMAX = ',nomax
1238	! ipause=mypauses(15,line)
1239	! call dadeb !(31,'ERR DANOT ',1)
1240	! endif
1241	!
1242	! not=nocut
1243	!
1244	! return
1245	! end subroutine getdanot
1246	subroutine daeps(deps)
1247	implicit none
1248	! **********************
1249	!
1250	! THIS SUBROUTINE RESETS THE TRUNCATION ORDER NOCUT TO A NEW VALUE
1251	!
1252	!-----------------------------------------------------------------------------
1253	!
1254	real(dp) deps
1255	!
1256	if(deps.ge.0.0_dp) then
1257	eps = deps
1258	else
1259	deps=eps
1260	endif
1261	!
1262	return
1263	end subroutine daeps
1264	!
1265	subroutine dapek(ina,jv,cjj)
1266	implicit none
1267	! ****************************
1268	!
1269	! THIS SUBROUTINE DETERMINES THE COEFFICIENT OF THE ARRAY
1270	! OF EXPONENTS JJ AND RETURNS IT IN CJJ
1271	!
1272	!-----------------------------------------------------------------------------
1273	!
1274	integer i,ibase,ic,ic1,ic2,icu,icz,ii_1,illa,ilma,ina,inoa,inva,ipek,ipoa,&
1275	iu,iz,jj1,mchk,ipause,mypauses
1276	integer,dimension(:)::jv ! 2002.12.4
1277	integer,dimension(lnv)::jj
1278	real(dp) cjj
1279	!
1280	if((.not.C_%STABLE_DA)) then
1281	if(c_%watch_user) then
1282	write(6,*) "big problem in dabnew ", sqrt(crash)
1283	endif
1284	return
1285	endif
1286	jj=0
1287	do i=1,size(jv)
1288	jj(i)=jv(i)
1289	enddo
1290	call dainf(ina,inoa,inva,ipoa,ilma,illa)
1291	if((.not.C_%STABLE_DA)) then
1292	if(c_%watch_user) then
1293	write(6,*) "big problem in dabnew ", sqrt(crash)
1294	endif
1295	return
1296	endif
1297	!
1298	!
1299	if(illa.eq.0) then ! Etienne shit
1300	cjj = 0
1301	return
1302	endif
1303	jj1 = 1
1304	if(inva.eq.0.or.nomax.eq.1) then
1305	if(inva.ne.0.and.nomax.eq.1) then
1306	if(illa.ge.2) then
1307	do i=1,illa - 1
1308	if(jj(i).eq.1) jj1 = i + 1
1309	enddo
1310	else
1311	jj1 = jj(1) + 1
1312	endif
1313	else
1314	jj1 = jj(1)
1315	endif
1316	if(jj1.lt.1.or.jj1.gt.illa) then
1317	print*,'ERROR IN DAPEK, INDEX OUTSIDE RANGE, JJ(1) = ',jj1
1318	! call dadeb !(31,'ERR DAPEK1',1)
1319	endif
1320	ipek = ipoa + jj1 - 1
1321	cjj = cc(ipek)
1322	return
1323	endif
1324	ii_1 = (nvmax+1)/2
1325	ibase = nomax+1
1326	!
1327	! DETERMINE INDEX TO BE SEARCHED FOR
1328	! **********************************
1329	!
1330	call dadcd(jj,ic1,ic2)
1331	!
1332	!ETIENNE
1333	if(ic1.gt.lia.or.ic2.gt.lia) then
1334	write(line,'(a24,i8)') 'DISASTER IN DAPEK, INA= ',ina
1335	ipause=mypauses(16,line)
1336	endif
1337	!ETIENNE
1338	ic = ia1(ic1) + ia2(ic2)
1339	!
1340	! DETERMINE IF MONOMIAL TO BE POKED CONFORMS WITH INOA, INVA,NOCUT
1341	! ****************************************************************
1342	!
1343	! IF(ICO.GT.INOA.OR.ICV.GT.INVA.OR.ICO.GT.NOCUT) THEN
1344	! CJJ = 0
1345	! RETURN
1346	! ENDIF
1347	!
1348	! DETERMINE IF MONOMIAL IS INSIDE FIRST AND LAST MONOMIALS OF A
1349	! *************************************************************
1350	!
1351	iu = ipoa
1352	iz = ipoa + illa - 1
1353	icu = ia1(i_1(iu))+ia2(i_2(iu))
1354	icz = ia1(i_1(iz))+ia2(i_2(iz))
1355	!
1356	if(illa.eq.0) then
1357	cjj = 0
1358	return
1359	elseif(ic.eq.icu) then
1360	cjj = cc(iu)
1361	return
1362	elseif(ic.eq.icz) then
1363	cjj = cc(iz)
1364	return
1365	elseif(ic.lt.icu.or.ic.gt.icz) then
1366	cjj = 0
1367	return
1368	endif
1369	!
1370	! SEARCHING PROPER MONOMIAL
1371	! *************************
1372	!
1373	10 continue
1374	if(iz-iu.le.1) then
1375	cjj = 0
1376	return
1377	endif
1378	i = (iu+iz)/2
1379	!
1380	! if(ia1(i_1(i))+ia2(i_2(i)) - ic) 20,30,40
1381	mchk=ia1(i_1(i))+ia2(i_2(i)) - ic
1382	if(mchk.lt.0) goto 20
1383	if(mchk.eq.0) goto 30
1384	if(mchk.gt.0) goto 40
1385	20 iu = i
1386	goto 10
1387	30 cjj = cc(i)
1388	return
1389	40 iz = i
1390	goto 10
1391	!
1392	end subroutine dapek
1393	!
1394	subroutine dapok(ina,jv,cjj)
1395	implicit none
1396	! ****************************
1397	!
1398	! THIS SUBROUTINE SETS THE COEFFICIENT OF THE ARRAY
1399	! OF EXPONENTS JJ TO THE VALUE CJJ
1400	!
1401	!-----------------------------------------------------------------------------
1402	!
1403	integer i,ic,ic1,ic2,icu,icz,ii,illa,ilma,ina,inoa,inva,ipoa,ipok,&
1404	iu,iz,jj1,mchk,ipause,mypauses
1405	integer,dimension(:)::jv ! 2002.12.4
1406	integer,dimension(lnv)::jj
1407	real(dp) cjj
1408	!
1409	if((.not.C_%STABLE_DA)) then
1410	if(c_%watch_user) then
1411	write(6,*) "big problem in dabnew ", sqrt(crash)
1412	endif
1413	return
1414	endif
1415	jj=0
1416	do i=1,size(jv)
1417	jj(i)=jv(i)
1418	enddo
1419	!
1420	call dainf(ina,inoa,inva,ipoa,ilma,illa)
1421	if((.not.C_%STABLE_DA)) then
1422	if(c_%watch_user) then
1423	write(6,*) "big problem in dabnew ", sqrt(crash)
1424	endif
1425	return
1426	endif
1427	!
1428	!
1429	jj1 = 1
1430	if(inva.eq.0.or.nomax.eq.1) then
1431	if(inva.ne.0.and.nomax.eq.1) then
1432	if(illa.ge.2) then
1433	do i=1,illa - 1
1434	if(jj(i).eq.1) jj1 = i + 1
1435	enddo
1436	else
1437	jj1 = jj(1) + 1
1438	endif
1439	else
1440	jj1 = jj(1)
1441	endif
1442	if(jj1.lt.1.or.jj1.gt.illa) then
1443	print*,'ERROR IN DAPOK, INDEX OUTSIDE RANGE, JJ(1) = ',jj1
1444	! call dadeb !(31,'ERR DAPOK1',1)
1445	endif
1446	ipok = ipoa + jj1 - 1
1447	cc(ipok) = cjj
1448	return
1449	endif
1450	! DETERMINE INDEX TO BE SEARCHED FOR
1451	! **********************************
1452	!
1453	call dadcd(jj,ic1,ic2)
1454	!
1455	ic = ia1(ic1) + ia2(ic2)
1456	!
1457	! DETERMINE IF MONOMIAL TO BE POKED CONFORMS WITH INOA, INVA,NOCUT
1458	! ****************************************************************
1459	!
1460	!
1461	if(illa.ne.0) then ! etienne shit
1462	iu = ipoa
1463	iz = ipoa + illa - 1
1464	!
1465	! DETERMINE IF MONOMIAL IS INSIDE FIRST AND LAST MONOMIALS OF A
1466	! *************************************************************
1467	!
1468	icu = ia1(i_1(iu))+ia2(i_2(iu))
1469	icz = ia1(i_1(iz))+ia2(i_2(iz))
1470	endif
1471	if(illa.eq.0) then
1472	i = ipoa
1473	goto 100
1474	elseif(ic.eq.icu) then
1475	cc(iu) = cjj
1476	i = iu
1477	goto 200
1478	elseif(ic.eq.icz) then
1479	cc(iz) = cjj
1480	i = iz
1481	goto 200
1482	elseif(ic.lt.icu) then
1483	i = iu
1484	goto 100
1485	elseif(ic.gt.icz) then
1486	i = iz + 1
1487	goto 100
1488	endif
1489	!
1490	!
1491	! SEARCHING PLACE TO POKE INTO OR BEFORE WHICH TO POKE
1492	! ****************************************************
1493	!
1494	iu = ipoa
1495	iz = ipoa + illa
1496	!
1497	10 continue
1498	if(iz-iu.le.1) then
1499	i = iz
1500	goto 100
1501	endif
1502	i = (iu+iz)/2
1503	!
1504	! if(ia1(i_1(i))+ia2(i_2(i)) - ic) 20,30,40
1505	mchk=ia1(i_1(i))+ia2(i_2(i)) - ic
1506	if(mchk.lt.0) goto 20
1507	if(mchk.eq.0) goto 30
1508	if(mchk.gt.0) goto 40
1509	20 iu = i
1510	goto 10
1511	30 cc(i) = cjj
1512	goto 200
1513	40 iz = i
1514	goto 10
1515	!
1516	! INSERTING THE MONOMIAL, MOVING THE REST
1517	! ***************************************
1518	!
1519	100 continue
1520	!
1521	if(abs(cjj).lt.eps) return
1522	!
1523	do ii=ipoa+illa,i+1,-1
1524	cc(ii) = cc(ii-1)
1525	i_1(ii) = i_1(ii-1)
1526	i_2(ii) = i_2(ii-1)
1527	enddo
1528	!
1529	cc(i) = cjj
1530	i_1(i) = ic1
1531	i_2(i) = ic2
1532	!
1533	idall(ina) = illa + 1
1534	if(idall(ina).gt.idalm(ina)) then
1535	write(line,'(a15)') 'ERROR IN DAPOK '
1536	ipause=mypauses(17,line)
1537	call dadeb !(31,'ERR DAPOK ',1)
1538	endif
1539	!
1540	return
1541	!
1542	! CASE OF CJJ = 0 WHICH MEANS MOVING THE REST
1543	! *********************************************
1544	!
1545	200 continue
1546	if(abs(cjj).lt.eps) then
1547	do ii=i,ipoa+illa-2
1548	cc(ii) = cc(ii+1)
1549	i_1(ii) = i_1(ii+1)
1550	i_2(ii) = i_2(ii+1)
1551	enddo
1552	idall(ina) = illa - 1
1553	endif
1554	return
1555	!
1556	end subroutine dapok
1557	!
1558	subroutine daclr(inc)
1559	implicit none
1560	! *********************
1561	!
1562	! THIS SUBROUTINE SETS ALL THE STACK SPACE RESERVED FOR VARIABLE
1563	! C TO ZERO
1564	!
1565	!-----------------------------------------------------------------------------
1566	!
1567	integer i,illc,ilmc,inc,inoc,invc,ipoc
1568	!
1569	if((.not.C_%STABLE_DA)) then
1570	if(c_%watch_user) then
1571	write(6,*) "big problem in dabnew ", sqrt(crash)
1572	endif
1573	return
1574	endif
1575	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
1576	if((.not.C_%STABLE_DA)) then
1577	if(c_%watch_user) then
1578	write(6,*) "big problem in dabnew ", sqrt(crash)
1579	endif
1580	return
1581	endif
1582	!
1583	do i=ipoc,ipoc+ilmc-1
1584	!
1585	cc(i) = 0.0_dp
1586	!
1587	enddo
1588	!
1589	return
1590	end subroutine daclr
1591	!
1592	subroutine dacop(ina,inb)
1593	implicit none
1594	! *************************
1595	!
1596	! THIS SUBROUTINE COPIES THE DA VECTOR A TO THE DA VECTOR B
1597	!
1598	!-----------------------------------------------------------------------------
1599	!
1600	! call dainf(ina,inoa,inva,ipoa,ilma,illa)
1601	! call dainf(inb,inob,invb,ipob,ilmb,illb)
1602	!
1603	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INB,INOB,INVB)
1604	!-----------------------------------------------------------------------------
1605	!
1606	integer ia,ib,illa,ina,inb,ipoa,ipob
1607	!
1608	if((.not.C_%STABLE_DA)) then
1609	if(c_%watch_user) then
1610	write(6,*) "big problem in dabnew ", sqrt(crash)
1611	endif
1612	return
1613	endif
1614	ipob = idapo(inb)
1615	ipoa = idapo(ina)
1616	illa = idall(ina)
1617	ib = ipob - 1
1618	!
1619	! iif = 0
1620	! if(nomax.eq.1.or.inva.eq.0) iif = 1
1621	do ia = ipoa,ipoa+illa-1
1622	!
1623	if(nomax.gt.1) then
1624	if(ieo(ia1(i_1(ia))+ia2(i_2(ia))).gt.nocut) goto 100
1625	endif
1626	ib = ib + 1
1627	cc(ib) = cc(ia)
1628	i_1(ib) = i_1(ia)
1629	i_2(ib) = i_2(ia)
1630	!
1631	100 continue
1632	enddo
1633	!
1634	idall(inb) = ib - ipob + 1
1635	return
1636	end subroutine dacop
1637	subroutine daadd(ina,inb,inc)
1638	implicit none
1639	! *****************************
1640	!
1641	! THIS SUBROUTINE PERFORMS A DA ADDITION OF THE DA VECTORS A AND B.
1642	! THE RESULT IS STORED IN C.
1643	!
1644	!-----------------------------------------------------------------------------
1645	!
1646	integer i,ina,ipoa
1647	integer idaadd,inb,inc,ipoc
1648	integer ipob
1649	!
1650	if((.not.C_%STABLE_DA)) then
1651	if(c_%watch_user) then
1652	write(6,*) "big problem in dabnew ", sqrt(crash)
1653	endif
1654	return
1655	endif
1656	if(nomax.eq.1) then
1657	ipoc = idapo(inc)
1658	ipoa = idapo(ina)
1659	ipob = idapo(inb)
1660	! minv = min(inva,invb,invc)
1661	do i=0,nvmax
1662	cc(ipoc+i) = cc(ipoa+i) + cc(ipob+i)
1663	enddo
1664	return
1665	endif
1666	if(ina.ne.inc.and.inb.ne.inc) then
1667	call dalin(ina,+1.0_dp,inb,+1.0_dp,inc)
1668	else
1669	idaadd = 0
1670	call daall1(idaadd,'$$DAADD $$',nomax,nvmax)
1671	call dalin(ina,+1.0_dp,inb,+1.0_dp,idaadd)
1672	call dacop(idaadd,inc)
1673	call dadal1(idaadd)
1674	endif
1675	!
1676	return
1677	end subroutine daadd
1678	!
1679	subroutine datrunc(ina,io,inb)
1680	implicit none
1681	integer ina,io,inb,nt
1682	if((.not.C_%STABLE_DA)) then
1683	if(c_%watch_user) then
1684	write(6,*) "big problem in dabnew ", sqrt(crash)
1685	endif
1686	return
1687	endif
1688	nt=nocut
1689	if(io>nomax) then
1690	if(ina/=inb) call dacop(ina,inb)
1691	return
1692	endif
1693	nocut=io-1
1694	call dacop(ina,inb)
1695	nocut = nt
1696	end subroutine datrunc
1697	subroutine dasub(ina,inb,inc)
1698	implicit none
1699	! THIS SUBROUTINE PERFORMS A DA SUBTRACTION OF THE DA VECTORS A AND B.
1700	! THE RESULT IS STORED IN C.
1701	!
1702	!-----------------------------------------------------------------------------
1703	!
1704	integer idasub
1705	integer i,ina,ipoa
1706	integer inc,ipoc,inb
1707	integer ipob
1708	!
1709	if((.not.C_%STABLE_DA)) then
1710	if(c_%watch_user) then
1711	write(6,*) "big problem in dabnew ", sqrt(crash)
1712	endif
1713	return
1714	endif
1715	if(nomax.eq.1) then
1716	ipoc = idapo(inc)
1717	ipoa = idapo(ina)
1718	ipob = idapo(inb)
1719	! minv = min(inva,invb,invc)
1720	do i=0,nvmax
1721	cc(ipoc+i) = cc(ipoa+i) - cc(ipob+i)
1722	enddo
1723	return
1724	endif
1725	if(ina.ne.inc.and.inb.ne.inc) then
1726	call dalin(ina,+1.0_dp,inb,-1.0_dp,inc)
1727	else
1728	idasub = -1
1729	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
1730	call daall1(idasub,'$$DASUB $$',nomax,nvmax)
1731	call dalin(ina,+1.0_dp,inb,-1.0_dp,idasub)
1732	call dacop(idasub,inc)
1733	call dadal1(idasub)
1734	endif
1735	!
1736	return
1737	end subroutine dasub
1738	subroutine damul(ina,inb,inc)
1739	implicit none
1740	! *****************************
1741	!
1742	! THIS SUBROUTINE PERFORMS A DA MULTIPLICATION OF THE DA VECTORS A AND B.
1743	! THE RESULT IS STORED IN C. AS TEMPORARY STORAGE, THE STACK SPACE
1744	! OF THE (NOMAX+2) SCRATCH VARIABLES ALLOCATED BY DAINI IS USED.
1745	!
1746	!-----------------------------------------------------------------------------
1747	!
1748	integer ina,inb,inc,incc,ipoc,ipoa,ipob,i
1749	real(dp) ccipoa,ccipob
1750	!
1751	if((.not.C_%STABLE_DA)) then
1752	if(c_%watch_user) then
1753	write(6,*) "big problem in dabnew ", sqrt(crash)
1754	endif
1755	return
1756	endif
1757	if(nomax.eq.1) then
1758	ipoa=idapo(ina)
1759	ipob=idapo(inb)
1760	ipoc=idapo(inc)
1761	! minv = min(inva,invb,invc)
1762	ccipoa = cc(ipoa)
1763	ccipob = cc(ipob)
1764	cc(ipoc) = ccipoa*ccipob
1765	do i=1,nvmax
1766	cc(ipoc+i) = ccipoacc(ipob+i) + ccipobcc(ipoa+i)
1767	enddo
1768	! do 30 i=ipoc+minv+1,ipoc+invc
1769	! 30 cc(i) = zero
1770	return
1771	endif
1772	if(ina.eq.inc.or.inb.eq.inc) then
1773	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
1774	incc=0
1775	call daall1(incc,'$$DAJUNK$$',nomax,nvmax)
1776	call damult(ina,inb,incc)
1777	call dacop(incc,inc)
1778	call dadal1(incc)
1779	else
1780	call damult(ina,inb,inc)
1781	endif
1782	return
1783	end subroutine damul
1784	subroutine damult(ina,inb,inc)
1785	implicit none
1786	! *****************************
1787	!
1788	! THIS SUBROUTINE PERFORMS A DA MULTIPLICATION OF THE DA VECTORS A AND B.
1789	! THE RESULT IS STORED IN C. AS TEMPORARY STORAGE, THE STACK SPACE
1790	! OF THE (NOMAX+2) SCRATCH VARIABLES ALLOCATED BY DAINI IS USED.
1791	!
1792	!-----------------------------------------------------------------------------
1793	!
1794	!
1795	!
1796	! CALL DACHK(INA,INOA,INVA, INB,INOB,INVB, INC,INOC,INVC)
1797	!
1798	! CASE OF FIRST ORDER ONLY
1799	! ************************
1800	!
1801	!-----------------------------------------------------------------------------
1802	!
1803	integer i,i_1ia,i_2ia,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,ina,inb,inc,inoa,inob,inoc,&
1804	inva,invb,invc,ioffb,ipoa,ipob,ipoc,ipos,noib,nom
1805	integer,dimension(0:lno)::ipno,noff
1806	real(dp) ccia,ccipoa,ccipob
1807	!
1808	if((.not.C_%STABLE_DA)) then
1809	if(c_%watch_user) then
1810	write(6,*) "big problem in dabnew ", sqrt(crash)
1811	endif
1812	return
1813	endif
1814	if(nomax.eq.1) then
1815	ipoa=idapo(ina)
1816	ipob=idapo(inb)
1817	ipoc=idapo(inc)
1818	! minv = min(inva,invb,invc)
1819	ccipoa = cc(ipoa)
1820	ccipob = cc(ipob)
1821	cc(ipoc) = ccipoa*ccipob
1822	do i=1,nvmax
1823	cc(ipoc+i) = ccipoacc(ipob+i) + ccipobcc(ipoa+i)
1824	enddo
1825	! do 30 i=ipoc+minv+1,ipoc+invc
1826	! 30 cc(i) = zero
1827	return
1828	endif
1829	call dainf(ina,inoa,inva,ipoa,ilma,illa)
1830	call dainf(inb,inob,invb,ipob,ilmb,illb)
1831	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
1832	if((.not.C_%STABLE_DA)) then
1833	if(c_%watch_user) then
1834	write(6,*) "big problem in dabnew ", sqrt(crash)
1835	endif
1836	return
1837	endif
1838	! GENERAL CASE
1839	! ************
1840	!
1841	do i=0,nomax
1842	noff(i) = idapo(i+2)
1843	ipno(i) = 0
1844	enddo
1845	!
1846	call daclr(1)
1847	!
1848	! RE-SORTING THE VECTOR B INTO PIECES THAT ARE OF ONLY ONE ORDER
1849	! *************************************************************
1850	!
1851	do ib=ipob,ipob+illb-1
1852	!
1853	noib = ieo(ia1(i_1(ib))+ia2(i_2(ib)))
1854	ipos = ipno(noib) + 1
1855	ipno(noib) = ipos
1856	inob = noff(noib) + ipos
1857	!
1858	cc(inob) = cc(ib)
1859	i_1(inob) = i_1(ib)
1860	i_2(inob) = i_2(ib)
1861	!
1862	enddo
1863	!
1864	do i=0,nomax
1865	idall(i+2) = ipno(i)
1866	enddo
1867	!
1868	! PERFORMING ACTUAL MULTIPLICATION
1869	! ********************************
1870	!
1871	nom = min(nocut,inoc)
1872	!
1873	do ia=ipoa,ipoa+illa-1
1874	!
1875	i_1ia = i_1(ia)
1876	i_2ia = i_2(ia)
1877	ccia = cc(ia)
1878	!
1879	do noib = 0,nom-ieo(ia1(i_1(ia))+ia2(i_2(ia)))
1880	!
1881	ioffb = noff(noib)
1882	!
1883	do ib = ioffb+1,ioffb+ipno(noib)
1884	!
1885	ic = ia2(i_2ia+i_2(ib)) + ia1(i_1ia + i_1(ib))
1886	! Georg says maybe needs if(ic/=0)
1887	if(ic/=0) then
1888	cc(ic) = cc(ic) + ccia*cc(ib)
1889	else
1890	write(6,*) " Georg warn me about ic could be 0.0_dp"
1891	stop 999
1892	endif
1893	!
1894	enddo
1895	enddo
1896	enddo
1897	!
1898	call dapac(inc)
1899	!
1900	return
1901	end subroutine damult
1902	!
1903	subroutine dadiv(ina,inb,inc)
1904	implicit none
1905	! *************************
1906	!
1907	! THIS SUBROUTINE SQUARES THE VECTOR A AND STORES THE RESULT IN C.
1908	!
1909	!-----------------------------------------------------------------------------
1910	!
1911	! THIS SUBROUTINE PERFORMS A DA DIVISION OF THE DA VECTORS A AND B.
1912	! THE RESULT IS STORED IN C.
1913	!
1914	!-----------------------------------------------------------------------------
1915	!
1916	integer idadiv,inb,ina,inc,ipoc,ipoa,ipob,i
1917	real(dp) ck,ck1
1918	!
1919	if((.not.C_%STABLE_DA)) then
1920	if(c_%watch_user) then
1921	write(6,*) "big problem in dabnew ", sqrt(crash)
1922	endif
1923	return
1924	endif
1925	if(nomax.eq.1) then
1926	! minv = min(inva,invb)
1927	ipoa = idapo(ina)
1928	ipob = idapo(inb)
1929	ipoc = idapo(inc)
1930	ck=1.0_dp/cc(ipob)
1931	ck1=cc(ipoa)*ck
1932	do i=1,nvmax
1933	cc(ipoc+i) = (cc(ipoa+i)-cc(ipob+i)ck1)ck
1934	enddo
1935	cc(ipoc)=ck1
1936	return
1937	endif
1938	idadiv = 0
1939	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
1940	call daall1(idadiv,'$$DADIV $$',nomax,nvmax)
1941	call dafun('INV ',inb,idadiv)
1942	call damul(ina,idadiv,inc)
1943	call dadal1(idadiv)
1944	!
1945	return
1946	end subroutine dadiv
1947	!
1948	subroutine dasqr(ina,inc)
1949	implicit none
1950	! *************************
1951	!
1952	! THIS SUBROUTINE SQUARES THE VECTOR A AND STORES THE RESULT IN C.
1953	!
1954	!-----------------------------------------------------------------------------
1955	!
1956	integer ina,inc,incc,ipoc,i,ipoa
1957	real(dp) ccipoa
1958	!
1959	! CASE OF FIRST ORDER ONLY
1960	! ************************
1961	if((.not.C_%STABLE_DA)) then
1962	if(c_%watch_user) then
1963	write(6,*) "big problem in dabnew ", sqrt(crash)
1964	endif
1965	return
1966	endif
1967	if(nomax.eq.1) then
1968	ipoc = idapo(inc)
1969	ipoa = idapo(ina)
1970	! minv = min(inva,invc)
1971	ccipoa = cc(ipoa)
1972	cc(ipoc) = ccipoa*ccipoa
1973	do i=1,nvmax
1974	cc(ipoc+i) = 2.0_dpccipoacc(ipoa+i)
1975	enddo
1976	return
1977	endif
1978	if(ina.eq.inc) then
1979	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
1980	incc=0
1981	call daall1(incc,'$$DAJUNK$$',nomax,nvmax)
1982	call dasqrt(ina,incc)
1983	call dacop(incc,inc)
1984	call dadal1(incc)
1985	else
1986	call dasqrt(ina,inc)
1987	endif
1988	return
1989	end subroutine dasqr
1990	subroutine dasqrt(ina,inc)
1991	implicit none
1992	! *************************
1993	!
1994	! THIS SUBROUTINE SQUARES THE VECTOR A AND STORES THE RESULT IN C.
1995	!
1996	!-----------------------------------------------------------------------------
1997	!
1998	!
1999	! CALL DACHK(INA,INOA,INVA,' ',-1,-1,INC,INOC,INVC)
2000	!
2001	!
2002	! CASE OF FIRST ORDER ONLY
2003	! ************************
2004	!
2005	!-----------------------------------------------------------------------------
2006	!
2007	integer i,i_1ia,i_2ia,ia,ib,ib1,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,&
2008	inva,invc,ioffa,ioffb,ipoa,ipoc,ipos,noia,noib,nom
2009	integer,dimension(0:lno)::ipno,noff
2010	real(dp) ccia,ccipoa
2011	!
2012	if((.not.C_%STABLE_DA)) then
2013	if(c_%watch_user) then
2014	write(6,*) "big problem in dabnew ", sqrt(crash)
2015	endif
2016	return
2017	endif
2018	if(nomax.eq.1) then
2019	ipoc = idapo(inc)
2020	ipoa = idapo(ina)
2021	! minv = min(inva,invc)
2022	ccipoa = cc(ipoa)
2023	cc(ipoc) = ccipoa*ccipoa
2024	do i=1,nvmax
2025	cc(ipoc+i) = 2.0_dpccipoacc(ipoa+i)
2026	enddo
2027	return
2028	endif
2029	call dainf(ina,inoa,inva,ipoa,ilma,illa)
2030	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
2031	if((.not.C_%STABLE_DA)) then
2032	if(c_%watch_user) then
2033	write(6,*) "big problem in dabnew ", sqrt(crash)
2034	endif
2035	return
2036	endif
2037	! GENERAL CASE
2038	! ************
2039	!
2040	do i=0,nomax
2041	noff(i) = idapo(i+2)
2042	ipno(i) = 0
2043	enddo
2044	!
2045	call daclr(1)
2046	!
2047	! RESORTING THE VECTOR A INTO PIECES THAT ARE OF ONLY ONE ORDER
2048	! *************************************************************
2049	!
2050	do ia=ipoa,ipoa+illa-1
2051	!
2052	noia = ieo(ia1(i_1(ia))+ia2(i_2(ia)))
2053	ipos = ipno(noia) + 1
2054	ipno(noia) = ipos
2055	inoa = noff(noia) + ipos
2056	!
2057	cc(inoa) = cc(ia)
2058	i_1(inoa) = i_1(ia)
2059	i_2(inoa) = i_2(ia)
2060	!
2061	enddo
2062	!
2063	do i=0,nomax
2064	idall(i+2) = ipno(i)
2065	enddo
2066	!
2067	! PERFORMING ACTUAL MULTIPLICATION
2068	! ********************************
2069	!
2070	nom = min(nocut,inoc)
2071	!
2072	do noia = 0,nom/2
2073	!
2074	ioffa = noff(noia)
2075	!
2076	do ia=ioffa+1,ioffa+ipno(noia)
2077	!
2078	i_1ia = i_1(ia)
2079	i_2ia = i_2(ia)
2080	ccia = cc(ia)
2081	!
2082	ic = ia2(i_2ia+i_2ia) + ia1(i_1ia+i_1ia)
2083	cc(ic) = cc(ic) + ccia*ccia
2084	ccia = ccia + ccia
2085	!
2086	do noib = noia,nom-noia
2087	!
2088	ioffb = noff(noib)
2089	if(noib.eq.noia) then
2090	ib1 = ia + 1
2091	else
2092	ib1 = ioffb + 1
2093	endif
2094	!
2095	do ib = ib1,ioffb+ipno(noib)
2096	!
2097	ic = ia2(i_2ia+i_2(ib)) + ia1(i_1ia + i_1(ib))
2098	cc(ic) = cc(ic) + ccia*cc(ib)
2099	!
2100	enddo
2101	enddo
2102	enddo
2103	enddo
2104	!
2105	call dapac(inc)
2106	!
2107	return
2108	end subroutine dasqrt
2109	!
2110	subroutine dacad(ina,ckon,inb)
2111	! use da_arrays
2112	implicit none
2113	! integer,dimension(lnv)::jjy
2114	! data jjy / lnv*0 / ! flat zero here
2115	! ******************************
2116	!
2117	! THIS SUBROUTINE ADDS THE CONSTANT CKON TO THE VECTOR A
2118	!
2119	!-----------------------------------------------------------------------------
2120	!
2121	integer ina,inb
2122	integer,parameter,dimension(lnv)::jjx=0
2123	real(dp) ckon,const
2124	!
2125	if((.not.C_%STABLE_DA)) then
2126	if(c_%watch_user) then
2127	write(6,*) "big problem in dabnew ", sqrt(crash)
2128	endif
2129	return
2130	endif
2131	call dacop(ina,inb)
2132	if(nomax.eq.1) then
2133	cc(idapo(inb)) = cc(idapo(inb)) + ckon
2134	return
2135	endif
2136	!
2137	call dapek(inb,jjx,const)
2138	call dapok(inb,jjx,const+ckon)
2139	!
2140	return
2141	end subroutine dacad
2142	!
2143	subroutine dacsu(ina,ckon,inb)
2144	implicit none
2145	! ******************************
2146	!
2147	! THIS SUBROUTINE SUBTRACTS THE CONSTANT CKON FROM THE VECTOR A
2148	!
2149	!-----------------------------------------------------------------------------
2150	!
2151	integer ina,inb
2152	integer,parameter,dimension(lnv)::jjx=0
2153	real(dp) ckon,const
2154	!
2155	if((.not.C_%STABLE_DA)) then
2156	if(c_%watch_user) then
2157	write(6,*) "big problem in dabnew ", sqrt(crash)
2158	endif
2159	return
2160	endif
2161	call dacop(ina,inb)
2162	!
2163	if(nomax.eq.1) then
2164	cc(idapo(inb)) = cc(idapo(inb)) - ckon
2165	return
2166	endif
2167	!
2168	call dapek(inb,jjx,const)
2169	call dapok(inb,jjx,const-ckon)
2170	!
2171	return
2172	end subroutine dacsu
2173	!
2174	subroutine dasuc(ina,ckon,inb)
2175	implicit none
2176	! ******************************
2177	!
2178	! THIS SUBROUTINE SUBTRACTS THE VECTOR INA FROM THE CONSTANT CKON
2179	!
2180	!-----------------------------------------------------------------------------
2181	!
2182	! call dainf(ina,inoa,inva,ipoa,ilma,illa)
2183	! call dainf(inb,inob,invb,ipob,ilmb,illb)
2184	!-----------------------------------------------------------------------------
2185	!
2186	integer i,ina,inb,ipoa,ipob
2187	real(dp) ckon
2188	!
2189	if((.not.C_%STABLE_DA)) then
2190	if(c_%watch_user) then
2191	write(6,*) "big problem in dabnew ", sqrt(crash)
2192	endif
2193	return
2194	endif
2195	ipob=idapo(inb)
2196	ipoa=idapo(ina)
2197	if(nomax.eq.1) then
2198	cc(ipob) = ckon - cc(ipoa)
2199	do i=1,nvmax
2200	cc(ipob+i) =-cc(ipoa+i)
2201	enddo
2202	return
2203	endif
2204	call dacsu(ina,ckon,inb)
2205	call dacmu(inb,-1.0_dp,inb)
2206	!
2207	return
2208	end subroutine dasuc
2209	!
2210	subroutine dacmu(ina,ckon,inc)
2211	implicit none
2212	! ******************************
2213	!
2214	! THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
2215	! THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
2216	! THE DA VECTOR DENOTED WITH THE INTEGER E.
2217	!
2218	!-----------------------------------------------------------------------------
2219	!
2220	integer ipoa,i,ina,inc,incc,ipoc
2221	real(dp) ckon
2222	!
2223	if((.not.C_%STABLE_DA)) then
2224	if(c_%watch_user) then
2225	write(6,*) "big problem in dabnew ", sqrt(crash)
2226	endif
2227	return
2228	endif
2229	if(nomax.eq.1) then
2230	! minv = min(inva,invb)
2231	ipoa = idapo(ina)
2232	ipoc = idapo(inc)
2233	do i=0,nvmax
2234	cc(ipoc+i) = cc(ipoa+i) * ckon
2235	enddo
2236	return
2237	endif
2238	if(ina.eq.inc) then
2239	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
2240	incc=0
2241	call daallno1(incc,'$$DAJUNK$$')
2242	call dacmut(ina,ckon,incc)
2243	call dacop(incc,inc)
2244	call dadal1(incc)
2245	else
2246	call dacmut(ina,ckon,inc)
2247	endif
2248	return
2249	end subroutine dacmu
2250	subroutine dacmut(ina,ckon,inb)
2251	implicit none
2252	! ******************************
2253	!
2254	! THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
2255	! THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
2256	! THE DA VECTOR DENOTED WITH THE INTEGER E.
2257	!
2258	!-----------------------------------------------------------------------------
2259	!
2260	!
2261	!
2262	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INB,INOB,INVB)
2263	!-----------------------------------------------------------------------------
2264	!
2265	integer i,ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,ipoa,&
2266	ipob,ipause,mypauses
2267	real(dp) ckon
2268	!
2269	if((.not.C_%STABLE_DA)) then
2270	if(c_%watch_user) then
2271	write(6,*) "big problem in dabnew ", sqrt(crash)
2272	endif
2273	return
2274	endif
2275	if(nomax.eq.1) then
2276	! minv = min(inva,invb)
2277	ipoa = idapo(ina)
2278	ipob = idapo(inb)
2279	do i=0,nvmax
2280	cc(ipob+i) = cc(ipoa+i) * ckon
2281	enddo
2282	return
2283	endif
2284	call dainf(ina,inoa,inva,ipoa,ilma,illa)
2285	call dainf(inb,inob,invb,ipob,ilmb,illb)
2286	if((.not.C_%STABLE_DA)) then
2287	if(c_%watch_user) then
2288	write(6,*) "big problem in dabnew ", sqrt(crash)
2289	endif
2290	return
2291	endif
2292	if(abs(ckon).lt.eps) then
2293	idall(inb) = 0
2294	return
2295	endif
2296	!
2297	ib = ipob - 1
2298	!
2299	do ia=ipoa,ipoa+illa-1
2300	!
2301	if(ieo(ia1(i_1(ia))+ia2(i_2(ia))).gt.nocut) goto 100
2302	ib = ib + 1
2303	cc(ib) = cc(ia)*ckon
2304	i_1(ib) = i_1(ia)
2305	i_2(ib) = i_2(ia)
2306	!
2307	100 continue
2308	enddo
2309	!
2310	idall(inb) = ib-ipob+1
2311	if(idall(inb).gt.idalm(inb)) then
2312	write(line,'(a17)') 'ERROR IN DACMU '
2313	ipause=mypauses(18,line)
2314	call dadeb !(31,'ERR DACMU ',1)
2315	endif
2316	!
2317	return
2318	end subroutine dacmut
2319	!
2320	subroutine dacdi(ina,ckon,inb)
2321	implicit none
2322	! ******************************
2323	!
2324	! THIS SUBROUTINE DIVIDES THE VECTOR INA BY THE CONSTANT CKON
2325	!
2326	!-----------------------------------------------------------------------------
2327	!
2328	integer i,ina,inb,ipoa,ipob,ipause,mypauses
2329	real(dp) ckon
2330	!
2331	if((.not.C_%STABLE_DA)) then
2332	if(c_%watch_user) then
2333	write(6,*) "big problem in dabnew ", sqrt(crash)
2334	endif
2335	return
2336	endif
2337	if(ckon==0.0_dp) then
2338	if(check_da) then
2339	C_%STABLE_DA=.false.
2340	messagelost='constant part 0.0_dp in dacdi'
2341	return
2342	else
2343	write(line,'(a38)') 'ERROR IN DACDI CKON IS 0.0_dp'
2344	ipause=mypauses(25,line)
2345	endif
2346	endif
2347	if(nomax.eq.1) then
2348	! minv = min(inva,invb)
2349	ipoa = idapo(ina)
2350	ipob = idapo(inb)
2351	do i=0,nvmax
2352	cc(ipob+i) = cc(ipoa+i)/ ckon
2353	enddo
2354	return
2355	endif
2356	!
2357	call dacmu(ina,1.0_dp/ckon,inb)
2358	!
2359	return
2360	end subroutine dacdi
2361	!
2362	!
2363	subroutine dadic(ina,ckon,inc)
2364	implicit none
2365	! ******************************
2366	!
2367	! THIS SUBROUTINE DIVIDES THE CONSTANT CKON BY THE VECTOR INA
2368	!
2369	!-----------------------------------------------------------------------------
2370	!
2371	integer i,idadic,ina,inc,ipoa,ipoc
2372	real(dp) ckon,ck
2373	!
2374	if((.not.C_%STABLE_DA)) then
2375	if(c_%watch_user) then
2376	write(6,*) "big problem in dabnew ", sqrt(crash)
2377	endif
2378	return
2379	endif
2380	ipoa = idapo(ina)
2381	if(cc(ipoa)==0.0_dp) then
2382	if(check_da) C_%STABLE_DA=.false.
2383	messagelost='constant part 0.0_dp in dadic'
2384	return
2385	endif
2386	if(nomax.eq.1) then
2387	! minv = min(inva,invb)
2388	ipoc = idapo(inc)
2389	cc(ipoc)=ckon/cc(ipoa)
2390	ck=cc(ipoc)/cc(ipoa)
2391	do i=1,nvmax
2392	cc(ipoc+i) = -cc(ipoa+i)* ck
2393	enddo
2394	return
2395	endif
2396	! if(abs(ckon).lt.eps) then !2002.11.28
2397	! call dacon(inc,zero)
2398	! return
2399	! endif
2400	idadic = 0
2401	call daall1(idadic,'$$DADIC $$',nomax,nvmax)
2402	! if(ckon.eq.zero) then
2403	! write(line,'(a18)') 'ERROR IN DACDI and DADIC, CKON IS ZERO' !2002.11.28
2404	! ipause=mypauses(19,line)
2405	! call dadeb !(31,'ERR DACDI ',1)
2406	! endif
2407	! call dacdi(ina,ckon,idadic)
2408	! call dafun('INV ',idadic,inc)
2409	call dafun('INV ',ina,idadic)
2410	call dacmu(idadic,ckon,inc)
2411	call dadal1(idadic)
2412	!
2413	return
2414	end subroutine dadic
2415	!
2416	subroutine dacma(ina,inb,bfac,inc)
2417	implicit none
2418	! **********************************
2419	!
2420	! THIS SUBROUTINE PERFORMS THE OPERATIONS C = A + B*BFAC, WHERE A,B,C ARE
2421	! DA VECTORS AND BFAC IS A real(dp). A AND C CAN BE IDENTICAL.
2422	! CAN LATER BE REPLACED BY SOMETHING LIKE DAADD WITH MINOR CHANGES.
2423	!
2424	!-----------------------------------------------------------------------------
2425	!
2426	integer idacma,ina,inb,inc,ipoc,ipob,ipoa,i
2427	real(dp) bfac
2428	!
2429	if((.not.C_%STABLE_DA)) then
2430	if(c_%watch_user) then
2431	write(6,*) "big problem in dabnew ", sqrt(crash)
2432	endif
2433	return
2434	endif
2435	if(nomax.eq.1) then
2436	ipoc = idapo(inc)
2437	ipoa = idapo(ina)
2438	ipob = idapo(inb)
2439	! minv = min(inva,invb,invc)
2440	do i=0,nvmax
2441	cc(ipoc+i) = cc(ipoa+i) + cc(ipob+i) * bfac
2442	enddo
2443	! do 8 i=ipoc+minv+1,ipoc+invc
2444	! 8 cc(i) = zero
2445	return
2446	endif
2447	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
2448	idacma = 0
2449	call daall1(idacma,'$$DACMA $$',nomax,nvmax)
2450	call dalin(ina,+1.0_dp,inb,bfac,idacma)
2451	call dacop(idacma,inc)
2452	call dadal1(idacma)
2453	!
2454	return
2455	end subroutine dacma
2456	!
2457	subroutine dalin(ina,afac,inb,bfac,inc)
2458	implicit none
2459	integer ina,inb,inc,incc,ipoc
2460	real(dp) afac,bfac
2461	! ***************************************
2462	!
2463	! THIS SUBROUTINE COMPUTES THE LINEAR COMBINATION
2464	! C = AFACA + BFACB. IT IS ALSO USED TO ADD AND SUBTRACT.
2465	!
2466	!-----------------------------------------------------------------------------
2467	!
2468	integer ipob,ipoa,i
2469	!
2470	if((.not.C_%STABLE_DA)) then
2471	if(c_%watch_user) then
2472	write(6,*) "big problem in dabnew ", sqrt(crash)
2473	endif
2474	return
2475	endif
2476	if(nomax.eq.1) then
2477	ipoc = idapo(inc)
2478	ipoa = idapo(ina)
2479	ipob = idapo(inb)
2480	! minv = min(inva,invb,invc)
2481	do i=0,nvmax
2482	cc(ipoc+i) = cc(ipoa+i) * afac + cc(ipob+i) * bfac
2483	enddo
2484	! do 8 i=ipoc+minv+1,ipoc+invc
2485	! 8 cc(i) = zero
2486	return
2487	endif
2488	if(ina.eq.inc.or.inb.eq.inc) then
2489	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
2490	incc=0
2491	call daall1(incc,'$$DAJUNK$$',nomax,nvmax)
2492	call dalint(ina,afac,inb,bfac,incc)
2493	call dacop(incc,inc)
2494	call dadal1(incc)
2495	else
2496	call dalint(ina,afac,inb,bfac,inc)
2497	endif
2498	return
2499	end subroutine dalin
2500	subroutine dalint(ina,afac,inb,bfac,inc)
2501	implicit none
2502	! ***************************************
2503	!
2504	! THIS SUBROUTINE COMPUTES THE LINEAR COMBINATION
2505	! C = AFACA + BFACB. IT IS ALSO USED TO ADD AND SUBTRACT.
2506	!
2507	!-----------------------------------------------------------------------------
2508	!
2509	! CALL DACHK(INA,INOA,INVA, INB,INOB,INVB, INC,INOC,INVC)
2510	!
2511	!-----------------------------------------------------------------------------
2512	!
2513	integer i,ia,iamax,ib,ibmax,ic,icmax,illa,illb,illc,ilma,ilmb,ilmc,ina,inb,inc,inoa,&
2514	inob,inoc,inva,invb,invc,ipoa,ipob,ipoc,is,ismax,ismin,ja,jb,mchk,&
2515	ipause,mypauses
2516	real(dp) afac,bfac,ccc,copf
2517	!
2518	if((.not.C_%STABLE_DA)) then
2519	if(c_%watch_user) then
2520	write(6,*) "big problem in dabnew ", sqrt(crash)
2521	endif
2522	return
2523	endif
2524	if(nomax.eq.1) then
2525	ipoc = idapo(inc)
2526	ipoa = idapo(ina)
2527	ipob = idapo(inb)
2528	! minv = min(inva,invb,invc)
2529	do i=0,nvmax
2530	cc(ipoc+i) = cc(ipoa+i) * afac + cc(ipob+i) * bfac
2531	enddo
2532	! do 8 i=ipoc+minv+1,ipoc+invc
2533	! 8 cc(i) = zero
2534	return
2535	endif
2536	call dainf(ina,inoa,inva,ipoa,ilma,illa)
2537	call dainf(inb,inob,invb,ipob,ilmb,illb)
2538	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
2539	if((.not.C_%STABLE_DA)) then
2540	if(c_%watch_user) then
2541	write(6,*) "big problem in dabnew ", sqrt(crash)
2542	endif
2543	return
2544	endif
2545	ia = ipoa
2546	ib = ipob
2547	ic = ipoc - 1
2548	iamax = ipoa+illa-1
2549	ibmax = ipob+illb-1
2550	icmax = ipoc+ilmc-1
2551	ja = ia1(i_1(ia)) + ia2(i_2(ia))
2552	jb = ia1(i_1(ib)) + ia2(i_2(ib))
2553	!
2554	if(ia.gt.iamax) then
2555	ismin = ib
2556	ismax = ibmax
2557	copf = bfac
2558	goto 50
2559	endif
2560	if(ib.gt.ibmax) then
2561	ismin = ia
2562	ismax = iamax
2563	copf = afac
2564	goto 50
2565	endif
2566	!
2567	! COMPARING
2568	! *********
2569	!
2570	10 continue
2571	! if(ja-jb) 30,20,40
2572	mchk=ja-jb
2573	if(mchk.lt.0) goto 30
2574	if(mchk.eq.0) goto 20
2575	if(mchk.gt.0) goto 40
2576	!
2577	! ADDING TWO TERMS
2578	! ****************
2579	!
2580	20 continue
2581	ccc = cc(ia)afac + cc(ib)bfac
2582	if(abs(ccc).lt.eps) goto 25
2583	if(ieo(ia1(i_1(ia))+ia2(i_2(ia))).gt.nocut) goto 25
2584	ic = ic + 1
2585	cc(ic) = ccc
2586	i_1(ic) = i_1(ia)
2587	i_2(ic) = i_2(ia)
2588	25 continue
2589	ia = ia + 1
2590	ib = ib + 1
2591	if(ia.gt.iamax) then
2592	ismin = ib
2593	ismax = ibmax
2594	copf = bfac
2595	goto 50
2596	endif
2597	if(ib.gt.ibmax) then
2598	ismin = ia
2599	ismax = iamax
2600	copf = afac
2601	goto 50
2602	endif
2603	ja = ia1(i_1(ia)) + ia2(i_2(ia))
2604	jb = ia1(i_1(ib)) + ia2(i_2(ib))
2605	goto 10
2606	!
2607	! STORING TERM A
2608	! **************
2609	!
2610	30 continue
2611	if(ieo(ia1(i_1(ia))+ia2(i_2(ia))).gt.nocut) goto 35
2612	ccc = cc(ia)*afac
2613	if(abs(ccc).lt.eps) goto 35
2614	ic = ic + 1
2615	cc(ic) = ccc
2616	i_1(ic) = i_1(ia)
2617	i_2(ic) = i_2(ia)
2618	35 continue
2619	ia = ia + 1
2620	if(ia.gt.iamax) then
2621	ismin = ib
2622	ismax = ibmax
2623	copf = bfac
2624	goto 50
2625	endif
2626	ja = ia1(i_1(ia)) + ia2(i_2(ia))
2627	goto 10
2628	!
2629	! STORING TERM B
2630	! **************
2631	!
2632	40 continue
2633	if(ieo(ia1(i_1(ib))+ia2(i_2(ib))).gt.nocut) goto 45
2634	ccc = cc(ib)*bfac
2635	if(abs(ccc).lt.eps) goto 45
2636	ic = ic + 1
2637	cc(ic) = ccc
2638	i_1(ic) = i_1(ib)
2639	i_2(ic) = i_2(ib)
2640	45 continue
2641	ib = ib + 1
2642	if(ib.gt.ibmax) then
2643	ismin = ia
2644	ismax = iamax
2645	copf = afac
2646	goto 50
2647	endif
2648	jb = ia1(i_1(ib)) + ia2(i_2(ib))
2649	goto 10
2650	!
2651	! COPYING THE REST
2652	! ****************
2653	!
2654	50 continue
2655	do is=ismin,ismax
2656	if(ieo(ia1(i_1(is))+ia2(i_2(is))).gt.nocut) goto 60
2657	ccc = cc(is)*copf
2658	if(abs(ccc).lt.eps) goto 60
2659	ic = ic + 1
2660	cc(ic) = ccc
2661	i_1(ic) = i_1(is)
2662	i_2(ic) = i_2(is)
2663	60 continue
2664	enddo
2665	!
2666	idall(inc) = ic - ipoc + 1
2667	!
2668	if(idall(inc).gt.idalm(inc)) then
2669	write(line,'(a40)') 'ERROR IN dalint idall(inc).gt.idalm(inc)'
2670	ipause=mypauses(21,line)
2671	call dadeb !(31,'ERR DALIN ',1)
2672	endif
2673	!
2674	return
2675	end subroutine dalint
2676	!
2677	subroutine dafun(cf,ina,inc)
2678	implicit none
2679	! ****************************
2680	!
2681	! THIS SUBROUTINE COMPUTES THE FUNCTION CF OF THE DA VECTOR A
2682	! AND STORES THE RESULT IN C.
2683	! AT PRESENT, SOME FUNCTIONS CAN BE COMPUTED ONLY TO FIFTH ORDER.
2684	! THIS HAS TO BE FIXED IN THE FUTURE.
2685	!
2686	!-----------------------------------------------------------------------------
2687	!
2688	integer ina,inc,incc
2689	character(4) cf
2690	!
2691	if((.not.C_%STABLE_DA)) then
2692	if(c_%watch_user) then
2693	write(6,*) "big problem in dabnew ", sqrt(crash)
2694	endif
2695	return
2696	endif
2697	if(ina.eq.inc) then
2698	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
2699	incc=0
2700	call daall1(incc,'$$DAJUNK$$',nomax,nvmax)
2701	call dafunt(cf,ina,incc)
2702	call dacop(incc,inc)
2703	call dadal1(incc)
2704	else
2705	call dafunt(cf,ina,inc)
2706	endif
2707	return
2708	end subroutine dafun
2709	subroutine dafunt(cf,ina,inc)
2710	implicit none
2711	! ****************************
2712	!
2713	! THIS SUBROUTINE COMPUTES THE FUNCTION CF OF THE DA VECTOR A
2714	! AND STORES THE RESULT IN C.
2715	! AT PRESENT, SOME FUNCTIONS CAN BE COMPUTED ONLY TO FIFTH ORDER.
2716	! THIS HAS TO BE FIXED IN THE FUTURE.
2717	!
2718	!-----------------------------------------------------------------------------
2719	!
2720	integer i,ina,inc,ind,inon,ipow,iscr,lfun,no,ipause,mypauses
2721	integer,parameter,dimension(lnv)::jjx=0
2722	real(dp) a0,ca,ea,ra,sa
2723	real(dp),dimension(0:lno)::xf
2724	character(4) cf,cfh
2725	character(26) abcs,abcc
2726	!
2727	data abcs /'abcdefghijklmnopqrstuvwxyz'/
2728	data abcc /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
2729	!
2730	if((.not.C_%STABLE_DA)) then
2731	if(c_%watch_user) then
2732	write(6,*) "big problem in dabnew ", sqrt(crash)
2733	endif
2734	return
2735	endif
2736	if(cf(1:1).eq.' ') then
2737	cfh(1:3) = cf(2:4)
2738	cfh(1:4) = ' '
2739	cf = cfh
2740	endif
2741	!
2742	do i=1,4
2743	ind = index(abcs,cf(i:i))
2744	if(ind.ne.0) cf(i:i) = abcc(ind:ind)
2745	enddo
2746	! call dainf(ina,inoa,inva,ipoa,ilma,illa)
2747	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
2748	!
2749	! CASE OF NV = 0 WHICH MEANS COORDINATEWISE OPERATION
2750	! ***************************************************
2751	!
2752	! CASE OF NV > 0 WHICH MEANS DIFFERENTIAL ALGEBRAIC OPERATION
2753	! ***********************************************************
2754	!
2755	if(cf.eq.'SQR ') then
2756	call dasqr(ina,inc)
2757	return
2758	endif
2759	!
2760	! ALLOCATE VARIABLES, PICK ZEROTH ORDER TERM
2761	! ******************************************
2762	!
2763	ipow = 0
2764	inon = 0
2765	iscr = 0
2766	!
2767	call daall1(ipow,'$$DAFUN1$$',nomax,nvmax)
2768	call daall1(inon,'$$DAFUN2$$',nomax,nvmax)
2769	call daall1(iscr,'$$DAFUN3$$',nomax,nvmax)
2770	!
2771	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC)
2772	!
2773	call dapek(ina,jjx,a0)
2774	!
2775	! no = min(nocut,inoa,inoc)
2776	no = min(nocut,nomax)
2777	!
2778	! BRANCHING TO DIFFERENT FUNCTIONS
2779	! ********************************
2780	!
2781	select case(cf)
2782	case('INV ')
2783	! if(cf.eq.'INV ') then
2784	! 1/(A0+P) = 1/A0(1-(P/A0)+(P/A0)*2-...)
2785	if(a0.eq.0) then
2786	if(check_da) then
2787	messagelost="a0.eq.0 for INV in dafun"
2788	C_%STABLE_DA=.false.
2789	C_%check_stable=.false.
2790	call dadal1(iscr)
2791	call dadal1(inon)
2792	call dadal1(ipow)
2793	return
2794	else
2795	write(*,1000) cf,ina,a0
2796	call dadeb !(31,'ERR DAFUN ',1)
2797	lfun = 0
2798	return
2799	endif
2800	endif
2801	xf(0) = 1.0_dp/a0
2802	do i=1,no
2803	xf(i) = -xf(i-1)/a0
2804	enddo
2805	!
2806	case('SQRT')
2807	! elseif(cf.eq.'SQRT') then
2808	! SQRT(A0+P) = SQRT(A0)(1+1/2(P/A0)-1/8(P/A0)**2+...)
2809	if(a0.le.0) then
2810	if(check_da) then
2811	messagelost="a0.le.0 for SQRT in dafun"
2812	C_%STABLE_DA=.false.
2813	C_%check_stable=.false.
2814	call dadal1(iscr)
2815	call dadal1(inon)
2816	call dadal1(ipow)
2817	return
2818	else
2819	write(*,1000) cf,ina,a0
2820	call dadeb !(31,'ERR DAFUN ',1)
2821	lfun = 0
2822	return
2823	endif
2824	endif
2825	ra = SQRT(a0)
2826	xf(0) = ra
2827	do i=1,no
2828	xf(i) = -xf(i-1)/a0/REAL(2i,kind=DP)REAL(2*i-3,kind=DP)
2829	enddo
2830	!
2831	!
2832	case('EXP ')
2833	! elseif(cf.eq.'EXP ') then
2834	! EXP(A0+P) = EXP(A0)(1+P+P*2/2!+...)
2835	if(a0>hyperbolic_aperture) then
2836	if(check_da) then
2837	messagelost="a0>hyperbolic_aperture for EXP in dafun"
2838	C_%STABLE_DA=.false.
2839	C_%check_stable=.false.
2840	call dadal1(iscr)
2841	call dadal1(inon)
2842	call dadal1(ipow)
2843	return
2844	else
2845	write(*,1000) cf,ina,a0
2846	call dadeb !(31,'ERR DAFUN ',1)
2847	lfun = 0
2848	return
2849	endif
2850	endif
2851	ea = exp(a0)
2852	xf(0) = ea
2853	do i=1,no
2854	xf(i) = xf(i-1)/REAL(i,kind=DP)
2855	enddo
2856	!
2857	case('LOG ')
2858	! elseif(cf.eq.'LOG ') then
2859	! LOG(A0+P) = LOG(A0) + (P/A0) - 1/2(P/A0)2 + 1/3(P/A0)**3 - ...)
2860	if(a0.le.0) then
2861	if(check_da) then
2862	messagelost="a0.le.0 for LOG in dafun"
2863	C_%STABLE_DA=.false.
2864	C_%check_stable=.false.
2865	call dadal1(iscr)
2866	call dadal1(inon)
2867	call dadal1(ipow)
2868	return
2869	else
2870	write(*,1000) cf,ina,a0
2871	call dadeb !(31,'ERR DAFUN ',1)
2872	lfun = 0
2873	return
2874	endif
2875	endif
2876	ea = LOG(a0)
2877	xf(0) = ea
2878	xf(1) = 1.0_dp/a0
2879	do i=2,no
2880	xf(i) = -xf(i-1)/a0/REAL(i,kind=DP)*REAL(i-1,kind=DP)
2881	enddo
2882	!
2883	case('SIN ')
2884	! elseif(cf.eq.'SIN ') then
2885	! SIN(A0+P) = SIN(A0)(1-P2/2!+P4/4!) + COS(A0)(P-P3/3!+P5/5!)
2886	sa = SIN(a0)
2887	ca = COS(a0)
2888	xf(0) = sa
2889	xf(1) = ca
2890	do i=2,no
2891	xf(i) = -xf(i-2)/REAL(i*(i-1),kind=DP)
2892	enddo
2893	!
2894	case('COS ')
2895	! elseif(cf.eq.'COS ') then
2896	! COS(A0+P) = COS(A0)(1-P2/2!+P4/4!) - SIN(A0)(P-P3/3!+P5/5!)
2897	sa = SIN(a0)
2898	ca = COS(a0)
2899	xf(0) = ca
2900	xf(1) = -sa
2901	do i=2,no
2902	xf(i) = -xf(i-2)/REAL(i*(i-1),kind=DP)
2903	enddo
2904	!
2905	case('SINH')
2906	! elseif(cf.eq.'SINH') then
2907	if(a0>hyperbolic_aperture) then
2908	if(check_da) then
2909	messagelost="a0>hyperbolic_aperture for SINH in dafun"
2910	C_%STABLE_DA=.false.
2911	C_%check_stable=.false.
2912	call dadal1(iscr)
2913	call dadal1(inon)
2914	call dadal1(ipow)
2915	return
2916	else
2917	write(*,1000) cf,ina,a0
2918	call dadeb !(31,'ERR DAFUN ',1)
2919	lfun = 0
2920	return
2921	endif
2922	endif
2923	sa = SINH(a0)
2924	ca = COSH(a0)
2925	xf(0) = sa
2926	xf(1) = ca
2927	do i=2,no
2928	xf(i) = xf(i-2)/REAL(i*(i-1),kind=DP)
2929	enddo
2930	case('COSH')
2931	! elseif(cf.eq.'COSH') then
2932	if(a0>hyperbolic_aperture) then
2933	if(check_da) then
2934	messagelost="a0>hyperbolic_aperture for COSH in dafun"
2935	C_%STABLE_DA=.false.
2936	C_%check_stable=.false.
2937	call dadal1(iscr)
2938	call dadal1(inon)
2939	call dadal1(ipow)
2940	return
2941	else
2942	write(*,1000) cf,ina,a0
2943	call dadeb !(31,'ERR DAFUN ',1)
2944	lfun = 0
2945	return
2946	endif
2947	endif
2948	sa = SINH(a0)
2949	ca = COSH(a0)
2950	xf(0) = ca
2951	xf(1) = sa
2952	xf(0) = ca
2953	xf(1) = sa
2954	do i=2,no
2955	xf(i) = xf(i-2)/REAL(i*(i-1),kind=DP)
2956	enddo
2957	case default
2958	! else
2959	write(line,'(a28,1x,a4)') 'ERROR, UNSOPPORTED FUNCTION ',cf
2960	ipause=mypauses(35,line)
2961	! endif
2962	end select
2963	!
2964	call dacon(inc,xf(0))
2965	call dacop(ina,inon)
2966	call dapok(inon,jjx,0.0_dp)
2967	call dacon(ipow,1.0_dp)
2968	!
2969	do i=1,min(no,nocut)
2970	!
2971	call damul(inon,ipow,iscr)
2972	call dacop(iscr,ipow)
2973	call dacma(inc,ipow,xf(i),inc)
2974	!
2975	enddo
2976	!
2977	1000 format('ERROR IN DAFUN, ',a4,' DOES NOT EXIST FOR VECTOR ',i10,'CONST TERM = ',e12.5)
2978	!
2979	call dadal1(iscr)
2980	call dadal1(inon)
2981	call dadal1(ipow)
2982	!
2983	return
2984	end subroutine dafunt
2985	!
2986	subroutine daabs(ina,anorm)
2987	implicit none
2988	! ***************************
2989	!
2990	! THIS SUBROUTINE COMPUTES THE NORM OF THE DA VECTOR A
2991	!
2992	!-----------------------------------------------------------------------------
2993	!
2994	integer i,illa,ilma,ina,inoa,inva,ipoa
2995	real(dp) anorm
2996	!
2997	if((.not.C_%STABLE_DA)) then
2998	if(c_%watch_user) then
2999	write(6,*) "big problem in dabnew ", sqrt(crash)
3000	endif
3001	return
3002	endif
3003	call dainf(ina,inoa,inva,ipoa,ilma,illa)
3004	if((.not.C_%STABLE_DA)) then
3005	if(c_%watch_user) then
3006	write(6,*) "big problem in dabnew ", sqrt(crash)
3007	endif
3008	return
3009	endif
3010	!
3011	anorm = 0.0_dp
3012	do i=ipoa,ipoa+illa-1
3013	anorm = anorm + abs(cc(i))
3014	enddo
3015	!
3016	return
3017	end subroutine daabs
3018	!
3019	!
3020	subroutine dacct(ma,ia,mb,ib,mc,ic)
3021	implicit none
3022	! ***********************************
3023	!
3024	! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC
3025	! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC
3026	! DA VECTORS, RESPECTIVELY.
3027	!
3028	!-----------------------------------------------------------------------------
3029	!
3030	integer i,ia,ib,ic,ij,illc,ilmc,inoc,invc,ipoc
3031	integer,dimension(lnv)::monx
3032	integer,dimension(:)::ma,mb,mc
3033	if((.not.C_%STABLE_DA)) then
3034	if(c_%watch_user) then
3035	write(6,*) "big problem in dabnew ", sqrt(crash)
3036	endif
3037	return
3038	endif
3039	!
3040	if(ma(1).eq.mc(1).or.mb(1).eq.mc(1)) then
3041	call dainf(mc(1),inoc,invc,ipoc,ilmc,illc)
3042	if((.not.C_%STABLE_DA)) then
3043	if(c_%watch_user) then
3044	write(6,*) "big problem in dabnew ", sqrt(crash)
3045	endif
3046	return
3047	endif
3048	do ij=1,ic
3049	monx(ij)=0
3050	enddo
3051	call daall(monx,ic,'$$DAJUNK$$',inoc,invc)
3052	call dacctt(ma,ia,mb,ib,monx,ic)
3053	do i=1,ic
3054	call dacop(monx(i),mc(i))
3055	enddo
3056	call dadal(monx,ic)
3057	else
3058	call dacctt(ma,ia,mb,ib,mc,ic)
3059	endif
3060	return
3061	end subroutine dacct
3062	subroutine dacctt(mb,ib,mc,ic,ma,ia)
3063	implicit none
3064	! ***********************************
3065	!
3066	! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC
3067	! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC
3068	! DA VECTORS, RESPECTIVELY.
3069	!
3070	!-----------------------------------------------------------------------------
3071	!
3072	! INTEGER MON(LNO+1),ICC(LNV)
3073	! INTEGER,dimension(:)::MB,MC,MA
3074	!ETIENNE
3075	!-----------------------------------------------------------------------------
3076	!
3077	integer i,ia,ib,ic,iia,iib,iic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc,&
3078	ipoa,ipob,ipoc,iv,jl,jv,ipause,mypauses
3079	integer,dimension(lno+1)::mon
3080	! integer,dimension(lno)::icc
3081	integer,dimension(lnv)::icc !johan 2008 March
3082	integer,dimension(:)::ma,mb,mc
3083	real(dp) ccf
3084	if((.not.C_%STABLE_DA)) then
3085	if(c_%watch_user) then
3086	write(6,*) "big problem in dabnew ", sqrt(crash)
3087	endif
3088	return
3089	endif
3090	!
3091	!ETIENNE
3092	!
3093	! CONSISTENCY CHECKS
3094	! ******************
3095	!
3096	iia = ma(1)
3097	iib = mb(1)
3098	iic = mc(1)
3099	call dainf(iia,inoa,inva,ipoa,ilma,illa)
3100	call dainf(iib,inob,invb,ipob,ilmb,illb)
3101	call dainf(iic,inoc,invc,ipoc,ilmc,illc)
3102	if((.not.C_%STABLE_DA)) then
3103	if(c_%watch_user) then
3104	write(6,*) "big problem in dabnew ", sqrt(crash)
3105	endif
3106	return
3107	endif
3108	!
3109	call damch(ma,ia)
3110	call damch(mb,ib)
3111	!
3112	if(ia.ne.ib) then
3113	write(line,'(a26)') 'ERROR IN DACCT, IA .NE. IB'
3114	ipause=mypauses(35,line)
3115	call dadeb !(31,'ERR DACCT1',1)
3116	elseif(ic.ne.invb) then
3117	write(line,'(a26)') 'ERROR IN DACCT, IC.NE.INVB'
3118	ipause=mypauses(35,line)
3119	call dadeb !(31,'ERR DACCT2',1)
3120	endif
3121	!
3122	! ALLOCATING LOCAL VECTORS AND CALLING MTREE
3123	! ******************************************
3124	!
3125	do i=1,ib
3126	icc(i) = 0
3127	enddo
3128	!
3129	do i=1,nomax+1
3130	mon(i) = 0
3131	enddo
3132	!
3133	call daall(icc,ib,'$$DACCT $$',nomax,nvmax)
3134	call daall(mon,nomax+1,'$$DAMON $$',inoc,invc)
3135	!
3136	call mtree(mb,ib,icc,ib)
3137	!
3138	! PERFORMING CONCATENATION
3139	! ************************
3140	!
3141	do i=1,ia
3142	call dacon(ma(i),cc(idapo(icc(i))))
3143	enddo
3144	!
3145	call dacon(mon(1),1.0_dp)
3146	!
3147	do i=1,idall(icc(1))-1
3148	!
3149	jl = i_1(idapo(icc(1))+i)
3150	jv = i_2(idapo(icc(1))+i)
3151	!
3152	call damul(mon(jl),mc(jv),mon(jl+1))
3153	!
3154	do iv=1,ia
3155	!
3156	ccf = cc(idapo(icc(iv))+i)
3157	if(abs(ccf).gt.eps) call dacma(ma(iv),mon(jl+1),ccf,ma(iv))
3158	!
3159	enddo
3160	enddo
3161	!
3162	call dadal(mon,nomax+1)
3163	call dadal(icc,ib)
3164	!
3165	return
3166	end subroutine dacctt
3167	!
3168	subroutine mtree(mb,ib,mc,ic)
3169	implicit none
3170	! *****************************
3171	!
3172	! THIS SUBROUTINE IS USED FOR CONCATENATION AND TRACKING OF VECTORS
3173	! THROUGH A DA MAP. IT COMPUTES THE TREE THAT HAS TO BE TRANSVERSED
3174	! MB IS THE DA MATRIX WITH IA TERMS. THE OUTPUT MC IS A CA MATRIX WHICH
3175	! CONTAINS COEFFICIENTS AND CONTROL INTEGERS USED FOR THE TRAVERSAL.
3176	!
3177	!-----------------------------------------------------------------------------
3178	!
3179	integer i,ib,ib1,ibi,ic,ic1,ic2,iccx,ichk,ii,iib,iic,illb,illc,ilmb,ilmc,&
3180	inob,inoc,invb,invc,ipob,ipoc,j,jl,jnon,nterm,ntermf,ipause,mypauses
3181	integer,dimension(lnv)::jjy
3182	integer,dimension(0:lno)::jv
3183	integer,dimension(:)::mb,mc
3184	real(dp) apek,bbijj,chkjj
3185	if((.not.C_%STABLE_DA)) then
3186	if(c_%watch_user) then
3187	write(6,*) "big problem in dabnew ", sqrt(crash)
3188	endif
3189	return
3190	endif
3191	!
3192	! CONSISTENCY CHECKS
3193	! ******************
3194	!
3195	iib = mb(1)
3196	iic = mc(1)
3197	call dainf(iib,inob,invb,ipob,ilmb,illb)
3198	call dainf(iic,inoc,invc,ipoc,ilmc,illc)
3199	if((.not.C_%STABLE_DA)) then
3200	if(c_%watch_user) then
3201	write(6,*) "big problem in dabnew ", sqrt(crash)
3202	endif
3203	return
3204	endif
3205	!
3206	call damch(mb,ib)
3207	call damch(mc,ic)
3208	!
3209	if(ib.ne.ic) then
3210	write(line,'(a26)') 'ERROR IN MTREE, IB .NE. IC'
3211	ipause=mypauses(35,line)
3212	call dadeb !(31,'ERR MTREE1',1)
3213	endif
3214	!
3215	! ALLOCATING LOCAL VECTORS
3216	! ************************
3217	!
3218	ichk = 0
3219	call daall1(ichk,'$$MTREE $$',nomax,nvmax)
3220	!
3221	! FIND ALL THE ENTRIES TO BE LOOKED FOR
3222	! *************************************
3223	!
3224	call daclr(1)
3225	!
3226	cc(1) = 1.0_dp
3227	!
3228	do i=1,ib
3229	if(nomax.eq.1) then
3230	do ib1 =2,invc+1 ! 2,7 ! Etienne
3231	cc(ib1) = 1.0_dp
3232	enddo
3233	else
3234	do ibi = idapo(mb(i)),idapo(mb(i))+idall(mb(i))-1
3235	iccx = ia1(i_1(ibi)) + ia2(i_2(ibi))
3236	if(ieo(iccx).gt.inob) goto 90
3237	cc(iccx) = 1.0_dp
3238	90 continue
3239	enddo
3240	endif
3241	enddo
3242	!
3243	do ii=1,inob
3244	!
3245	! SEARCHING FOR FATHER FOR EACH TERM
3246	!
3247	do i=1,nmmax
3248	if(cc(i).lt.0.5_dp) goto 140
3249	!
3250	jnon = 0
3251	call dancd(i_1(i),i_2(i),jjy)
3252	do j=1,invb
3253	if(jjy(j).eq.0) goto 130
3254	jnon = j
3255	jjy(j) = jjy(j) - 1
3256	call dadcd(jjy,ic1,ic2)
3257	apek = cc(ia1(ic1)+ia2(ic2))
3258	jjy(j) = jjy(j) + 1
3259	if(apek.ge.0.5_dp) goto 140
3260	130 continue
3261	enddo
3262	!
3263	if(jnon.eq.0) goto 140
3264	!
3265	! TERM IS AN ORPHAN, SO CREATE FOSTER FATHER
3266	!
3267	jjy(jnon) = jjy(jnon) - 1
3268	call dadcd(jjy,ic1,ic2)
3269	cc(ia1(ic1)+ia2(ic2)) = 1.0_dp
3270	!
3271	140 continue
3272	enddo
3273	enddo
3274	!
3275	call dapac(ichk)
3276	!ETIENNE CALL DAPRI(ICHK,32)
3277	!
3278	! SETTING UP TREE STRUCTURE
3279	! *************************
3280	!
3281	ntermf = idall(ichk)
3282	!
3283	! ZEROTH ORDER TERMS
3284	! ******************
3285	!
3286	do i=1,lnv
3287	jjy(i) = 0
3288	enddo
3289	!
3290	do i=1,ib
3291	call dapek(mb(i),jjy,bbijj)
3292	i_1(idapo(mc(i))) = 0
3293	i_2(idapo(mc(i))) = 0
3294	cc(idapo(mc(i))) = bbijj
3295	enddo
3296	!
3297	call dapek(ichk,jjy,chkjj)
3298	if(chkjj.gt.0.5_dp) then
3299	call dapok(ichk,jjy,-1.0_dp)
3300	else
3301	write(line,'(a49)') 'ERROR IN MTREE, ZEROTH ORDER TERM OF ICHK IS 0.0_dp'
3302	ipause=mypauses(35,line)
3303	call dadeb !(31,'ERR MTREE2',1)
3304	endif
3305	!
3306	nterm = 1
3307	!
3308	! HIGHER ORDER TERMS
3309	! ******************
3310	!
3311	do jl=1,inob
3312	jv(jl) = 0
3313	enddo
3314	!
3315	jl = 0
3316	chkjj = 1.0_dp
3317	!
3318	200 continue
3319	if(jl.eq.0.and.chkjj.le.0.5_dp) goto 250
3320	if(jl.lt.inob.and.chkjj.gt.0.5_dp) then
3321	jl = jl + 1
3322	jjy(1) = jjy(1) + 1
3323	jv(jl) = 1
3324	elseif(jv(jl).eq.invb) then
3325	jjy(jv(jl)) = jjy(jv(jl)) - 1
3326	jv(jl) = 0
3327	jl = jl - 1
3328	chkjj = 0.0_dp
3329	goto 200
3330	else
3331	jjy(jv(jl)) = jjy(jv(jl)) - 1
3332	jv(jl) = jv(jl) + 1
3333	jjy(jv(jl)) = jjy(jv(jl)) + 1
3334	endif
3335	!
3336	call dapek(ichk,jjy,chkjj)
3337	!
3338	if(chkjj.le.0.5_dp) goto 200
3339	!
3340	nterm = nterm + 1
3341	if(nterm.gt.idalm(mc(1))) then
3342	write(line,'(a31)') 'ERROR IN MTREE, NTERM TOO LARGE'
3343	ipause=mypauses(35,line)
3344	call dadeb !(31,'ERR MTREE3',1)
3345	endif
3346	!
3347	call dapok(ichk,jjy,-1.0_dp)
3348	!
3349	do i=1,ib
3350	call dapek(mb(i),jjy,bbijj)
3351	i_1(idapo(mc(i))+nterm-1) = jl
3352	i_2(idapo(mc(i))+nterm-1) = jv(jl)
3353	cc(idapo(mc(i))+nterm-1) = bbijj
3354	enddo
3355	!
3356	goto 200
3357	!
3358	250 continue
3359	!
3360	do i=1,ib
3361	idall(mc(i)) = nterm
3362	enddo
3363	!
3364	! PERFORMING CROSS CHECKS
3365	! ***********************
3366	!
3367	if(nterm.ne.ntermf.or.nterm.ne.idall(ichk)) then
3368	write(line,'(a46,1x,i8,1x,i8,1x,i8)') 'ERROR IN MTREE, NTERM, NTERMF, IDALL(ICHK) = ',nterm,ntermf,idall(ichk)
3369	ipause=mypauses(35,line)
3370	call dadeb !(31,'ERR MTREE4',1)
3371	endif
3372	!
3373	do i=idapo(ichk),idapo(ichk)+nterm-1
3374	if(abs(cc(i)+1.0_dp).gt.epsmac) then
3375	write(line,'(a44)') 'ERROR IN MTREE, NOT ALL TERMS IN ICHK ARE -1'
3376	ipause=mypauses(35,line)
3377	call dadeb !(31,'ERR MTREE5',1)
3378	endif
3379	enddo
3380	!
3381	call dadal1(ichk)
3382	!
3383	return
3384	end subroutine mtree
3385	!
3386	subroutine ppushprint(mc,ic,mf,jc,line)
3387	implicit none
3388	!
3389	integer i,ic,iv,jc,jl,jv,mf
3390	integer,dimension(:)::mc
3391	character(20) line
3392	if((.not.C_%STABLE_DA)) then
3393	if(c_%watch_user) then
3394	write(6,*) "big problem in dabnew ", sqrt(crash)
3395	endif
3396	return
3397	endif
3398	!
3399	if(mf.le.0) return
3400	write(mf,*) 0,0,jc+1,0,line
3401	do i=1,ic
3402	jc=1+jc
3403	write(mf,*) jc,jl,jv,cc(idapo(mc(i)))
3404	enddo
3405	! xf(i) = cc(idapo(mc(i)))
3406	! xm(1) = one
3407	do i=1,idall(mc(1))-1
3408	jl = i_1(idapo(mc(1))+i)
3409	jv = i_2(idapo(mc(1))+i)
3410	! xx = xm(jl)*xi(jv)
3411	! xm(jl+1) = xx
3412	do iv=1,ic
3413	jc=1+jc
3414	write(mf,*) jc,jl,jv,cc(idapo(mc(iv))+i)
3415	! xf(iv) = xf(iv) + cc(idapo(mc(iv))+i) * xx
3416	enddo
3417	enddo
3418	return
3419	end subroutine ppushprint
3420	!
3421	subroutine ppushstore(mc,nd2,coef,ml,mv)
3422	implicit none
3423	!
3424	integer i,ic,iv,jc,jl,jv,ntot,nd2
3425	integer,dimension(:), intent(in)::mc
3426	integer,dimension(:), intent(out)::ml,mv
3427	real(dp),dimension(:),intent(out)::coef
3428	!
3429	jc=0
3430	jl=0;jv=0;
3431	ic=nd2
3432	ntot=idall(mc(1))*ic
3433	do i=1,ic
3434	jc=1+jc
3435	ml(jc)=0
3436	mv(jc)=0
3437	coef(jc)=cc(idapo(mc(i)))
3438	! write(mf,*) jc,jl,jv,cc(idapo(mc(i)))
3439	enddo
3440	do i=1,idall(mc(1))-1
3441	jl = i_1(idapo(mc(1))+i)
3442	jv = i_2(idapo(mc(1))+i)
3443	! xx = xm(jl)*xi(jv)
3444	! xm(jl+1) = xx
3445	do iv=1,ic
3446	jc=1+jc
3447	ml(jc)=jl
3448	mv(jc)=jv
3449	coef(jc)=cc(idapo(mc(iv))+i)
3450	! write(mf,*) jc,jl,jv,cc(idapo(mc(iv))+i)
3451	! xf(iv) = xf(iv) + cc(idapo(mc(iv))+i) * xx
3452	enddo
3453	enddo
3454	return
3455	end subroutine ppushstore
3456	subroutine ppushGETN(mc,ND2,ntot)
3457	implicit none
3458	!
3459	integer ntot,ND2
3460	integer,dimension(:), intent(inout)::mc
3461	!
3462	ntot=idall(mc(1))*nd2
3463	end subroutine ppushGETN
3464	!
3465	subroutine ppush(mc,ic,xi,xf)
3466	implicit none
3467	! *****************************
3468	!
3469	! THIS SUBROUTINE APPLIES THE MATRIX WHOSE TREE IS STORED IN CA VECTOR MC
3470	! TO THE COORDINATES IN XI AND STORES THE RESULT IN XF
3471	!
3472	!-----------------------------------------------------------------------------
3473	!
3474	integer i,ic,iv,jl,jv
3475	integer,dimension(:)::mc
3476	real(dp) xx
3477	real(dp),dimension(lno+1)::xm
3478	real(dp),dimension(lno)::xt
3479	real(dp),dimension(:)::xf,xi
3480	if((.not.C_%STABLE_DA)) then
3481	if(c_%watch_user) then
3482	write(6,*) "big problem in dabnew ", sqrt(crash)
3483	endif
3484	return
3485	endif
3486	!
3487	do i=1,ic
3488	xt(i)=xi(i)
3489	enddo
3490	do i=1,ic
3491	xf(i) = cc(idapo(mc(i)))
3492	enddo
3493	!
3494	xm(1) = 1.0_dp
3495	!
3496	do i=1,idall(mc(1))-1
3497	!
3498	jl = i_1(idapo(mc(1))+i)
3499	jv = i_2(idapo(mc(1))+i)
3500	xx = xm(jl)*xt(jv)
3501	xm(jl+1) = xx
3502	!
3503	do iv=1,ic
3504	xf(iv) = xf(iv) + cc(idapo(mc(iv))+i) * xx
3505	enddo
3506	enddo
3507	do i=1,nvmax
3508	if(abs(xf(i))>da_absolute_aperture.and.check_da) then
3509	C_%STABLE_DA=.false.
3510	write(6,*) "unstable in ppush ",i,xf(i)
3511	endif
3512	enddo
3513	!
3514	return
3515	end subroutine ppush
3516	subroutine ppush1(mc,xi,xf)
3517	implicit none
3518	! *****************************
3519	!
3520	! THIS SUBROUTINE APPLIES THE MATRIX WHOSE TREE IS STORED IN CA VECTOR MC
3521	! TO THE COORDINATES IN XI AND STORES THE RESULT IN XF
3522	!
3523	!-----------------------------------------------------------------------------
3524	!
3525	integer i,jl,jv,mc
3526	real(dp) xf,xx
3527	real(dp),dimension(:)::xi
3528	real(dp),dimension(lno)::xt
3529	real(dp),dimension(lno+1)::xm
3530	if((.not.C_%STABLE_DA)) then
3531	if(c_%watch_user) then
3532	write(6,*) "big problem in dabnew ", sqrt(crash)
3533	endif
3534	return
3535	endif
3536	!
3537	do i=1,nvmax
3538	xt(i)=xi(i)
3539	enddo
3540	xf = cc(idapo(mc))
3541	!
3542	xm(1) = 1.0_dp
3543	!
3544	do i=1,idall(mc)-1
3545	!
3546	jl = i_1(idapo(mc)+i)
3547	jv = i_2(idapo(mc)+i)
3548	xx = xm(jl)*xt(jv)
3549	xm(jl+1) = xx
3550	!
3551	xf = xf + cc(idapo(mc)+i) * xx
3552	enddo
3553	if(abs(xf)>da_absolute_aperture.and.check_da) then
3554	C_%STABLE_DA=.false.
3555	write(6,*) "unstable in ppush1 ", xf
3556	write(6,*) xi(1:nvmax)
3557	endif
3558	return
3559	end subroutine ppush1
3560	subroutine dainv(ma,ia,mb,ib)
3561	implicit none
3562	! *****************************
3563	!
3564	! THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND
3565	! STORES THE RESULT IN MI
3566	!
3567	!-----------------------------------------------------------------------------
3568	!
3569	integer i,ia,ib,ij,illb,ilmb,inob,invb,ipob
3570	integer,dimension(lnv)::jj,ml
3571	integer,dimension(:)::ma,mb
3572	real(dp),dimension(lnv)::x
3573	if((.not.C_%STABLE_DA)) then
3574	if(c_%watch_user) then
3575	write(6,*) "big problem in dabnew ", sqrt(crash)
3576	endif
3577	return
3578	endif
3579	!
3580	do i=1,lnv
3581	jj(i)=0
3582	enddo
3583	if(ma(1).eq.mb(1)) then
3584	call dainf(mb(1),inob,invb,ipob,ilmb,illb)
3585	if((.not.C_%STABLE_DA)) then
3586	if(c_%watch_user) then
3587	write(6,*) "big problem in dabnew ", sqrt(crash)
3588	endif
3589	return
3590	endif
3591	do i=1,ia
3592	call dapok(ma(i),jj,0.0_dp)
3593	enddo
3594	do ij=1,ib
3595	ml(ij)=0
3596	enddo
3597	call daall(ml,ib,'$$DAJUNK$$',inob,invb)
3598	call dainvt(ma,ia,ml,ib)
3599	do i=1,ib
3600	call dacop(ml(i),mb(i))
3601	enddo
3602	call dadal(ml,ib)
3603	else
3604	do i=1,ia
3605	call dapek(ma(i),jj,x(i))
3606	call dapok(ma(i),jj,0.0_dp)
3607	enddo
3608	call dainvt(ma,ia,mb,ib)
3609	do i=1,ia
3610	call dapok(ma(i),jj,x(i))
3611	enddo
3612	endif
3613	return
3614	end subroutine dainv
3615	subroutine dainvt(ma,ia,mb,ib)
3616	implicit none
3617	! *****************************
3618	!
3619	! THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND
3620	! STORES THE RESULT IN MI
3621	!
3622	!-----------------------------------------------------------------------------
3623	!
3624	integer i,ia,ib,ie,ier,illa,illb,ilma,ilmb,inoa,inob,inva,invb,&
3625	ipoa,ipob,j,k,nocut0,ipause,mypauses
3626	integer,dimension(lnv)::jj,ms,ml
3627	integer,dimension(:)::ma,mb
3628	real(dp),dimension(lnv,lnv)::aa,ai
3629	real(dp) amjj,amsjj,prod
3630	if((.not.C_%STABLE_DA)) then
3631	if(c_%watch_user) then
3632	write(6,*) "big problem in dabnew ", sqrt(crash)
3633	endif
3634	return
3635	endif
3636	!
3637	call dainf(ma(1),inoa,inva,ipoa,ilma,illa)
3638	call dainf(mb(1),inob,invb,ipob,ilmb,illb)
3639	if((.not.C_%STABLE_DA)) then
3640	if(c_%watch_user) then
3641	write(6,*) "big problem in dabnew ", sqrt(crash)
3642	endif
3643	return
3644	endif
3645	!
3646	! CONSISTENCY CHECKS
3647	! ******************
3648	!
3649	call damch(ma,ia)
3650	call damch(mb,ib)
3651	!etienne
3652	do ie=1,ib
3653	call dacon(mb(ie),0.0_dp)
3654	enddo
3655	!etienne
3656	!
3657	if(ia.ne.ib) then
3658	write(line,'(a26)') 'ERROR IN DAINV, IA .NE. IB'
3659	ipause=mypauses(35,line)
3660	call dadeb !(31,'ERR DAINV1',1)
3661	elseif(ia.ne.inva.or.ib.ne.invb) then
3662	write(line,'(a40)') 'ERROR IN DAINV, IA.NE.INVA.OR.IB.NE.INVB'
3663	ipause=mypauses(35,line)
3664	call dadeb !(31,'ERR DAINV2',1)
3665	endif
3666	!
3667	! ALLOCATING LOCAL VECTORS
3668	! ************************
3669	!
3670	do i=1,ia
3671	ms(i) = 0
3672	ml(i) = 0
3673	enddo
3674	!
3675	call daall(ms,ia,'$$INV $$',inoa,inva)
3676	call daall(ml,ia,'$$INVL $$',inoa,inva)
3677	!
3678	! EXTRACTING LINEAR MATRIX, GENERATING NONLINEAR PART OF A
3679	! ********************************************************
3680	!
3681	do i=1,ib
3682	do j=1,ib
3683	do k=1,ib
3684	jj(k) = 0
3685	enddo
3686	jj(j) = 1
3687	call dapek(ma(i),jj,amjj)
3688	if(abs(amjj).gt.eps) call dapok(ma(i),jj,0.0_dp)
3689	aa(i,j) = amjj
3690	enddo
3691	call dacmu(ma(i),-1.0_dp,ma(i))
3692	enddo
3693	!
3694	! INVERTING LINEAR MATRIX, CHECKING RESULT AND STORING IN ML
3695	! **********************************************************
3696	!
3697	call matinv(aa,ai,ia,lnv,ier)
3698	!
3699	if(ier.eq.132) then
3700	if(check_da) then
3701	C_%STABLE_DA=.false.
3702	C_%check_stable=.false.
3703	messagelost='ERROR IN ROUTINE DAINV, ier=132 in matinv'
3704	write(6,*) messagelost
3705	call dadal(ml,ia)
3706	call dadal(ms,ia)
3707	return
3708	else
3709	write(line,'(a22)') 'ERROR IN ROUTINE DAINV'
3710	ipause=mypauses(35,line)
3711	call dadeb !(31,'ERR DAINV3',1)
3712	endif
3713	endif
3714	!
3715	ier = 0
3716	do i=1,ib
3717	do j=1,ib
3718	prod = 0.0_dp
3719	do k=1,ib
3720	jj(k) = 0
3721	prod = prod + aa(i,k)*ai(k,j)
3722	enddo
3723	if(i.eq.j) prod = prod - 1.0_dp
3724	if(abs(prod).gt.100.0_dp*epsmac) then
3725	write(6,) " abs(prod) > 100.0_dpepsmac in dainvt",abs(prod), 100.0_dp*epsmac
3726	if(check_da) then
3727	C_%STABLE_DA=.false.
3728	messagelost='ERROR IN ROUTINE DAINV, abs(prod).gt.100.0_dp*epsmac '
3729	call dadal(ml,ia)
3730	call dadal(ms,ia)
3731	return
3732	else
3733	write(line,'(a50,2(1x,i4),3(1x,g13.6))') 'ERROR IN DAINV, INVERSION DID NOT WORK,I,J,PROD = ' &
3734	& ,i,j,prod,epsmac,eps
3735	ipause=mypauses(35,line)
3736	ier = 1
3737	!ETIENNE
3738	return
3739	!ETIENNE
3740	endif
3741	endif
3742	jj(j) = 1
3743	call dapok(mb(i),jj,ai(i,j))
3744	call dapok(ml(i),jj,ai(i,j))
3745	enddo
3746	enddo
3747	!
3748	if(ier.eq.1) call dadeb !(31,'ERR DAINV4',1)
3749	!
3750	! ITERATIVELY COMPUTING DIFFERENT PARTS OF THE INVERSE
3751	! ****************************************************
3752	!
3753	! MB (OF ORDER I) = A1^-1 o [ E - ANL (NONLINEAR) o MB (OF ORDER I) ]
3754	!
3755	nocut0 = nocut
3756	!
3757	do i=2,nocut
3758	!
3759	nocut = i
3760	!
3761	call dacct(ma,ia,mb,ib,ms,ia)
3762	do j=1,ib
3763	do k=1,ib
3764	jj(k) = 0
3765	enddo
3766	jj(j) = 1
3767	call dapek(ms(j),jj,amsjj)
3768	call dapok(ms(j),jj,amsjj+1.0_dp)
3769	enddo
3770	!
3771	call dacct(ml,ia,ms,ia,mb,ib)
3772	!
3773	enddo
3774	!
3775	nocut = nocut0
3776	!
3777	! FLIPPING BACK SIGN OF A, FILLING UP FIRST ORDER PART AGAIN
3778	! **********************************************************
3779	!
3780	do i=1,ib
3781	call dacmu(ma(i),-1.0_dp,ma(i))
3782	do j=1,ib
3783	do k=1,ib
3784	jj(k) = 0
3785	enddo
3786	jj(j) = 1
3787	call dapok(ma(i),jj,aa(i,j))
3788	enddo
3789	enddo
3790	!
3791	call dadal(ml,ia)
3792	call dadal(ms,ia)
3793	!
3794	return
3795	end subroutine dainvt
3796	!
3797	subroutine dapin(ma,ia,mb,ib,jx)
3798	implicit none
3799	! *****************************
3800	!
3801	! THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND
3802	! STORES THE RESULT IN MI
3803	!
3804	!-----------------------------------------------------------------------------
3805	!
3806	integer i,ia,ib,ij,illb,ilmb,inob,invb,ipob
3807	integer,dimension(lnv)::jj,ml
3808	integer,dimension(:)::ma,mb,jx
3809	real(dp),dimension(lnv)::x
3810	!
3811	if((.not.C_%STABLE_DA)) then
3812	if(c_%watch_user) then
3813	write(6,*) "big problem in dabnew ", sqrt(crash)
3814	endif
3815	return
3816	endif
3817	do i=1,lnv
3818	jj(i)=0
3819	enddo
3820	if(ma(1).eq.mb(1)) then
3821	call dainf(mb(1),inob,invb,ipob,ilmb,illb)
3822	if((.not.C_%STABLE_DA)) then
3823	if(c_%watch_user) then
3824	write(6,*) "big problem in dabnew ", sqrt(crash)
3825	endif
3826	return
3827	endif
3828	do i=1,ia
3829	call dapok(ma(i),jj,0.0_dp)
3830	enddo
3831	do ij=1,ib
3832	ml(ij)=0
3833	enddo
3834	call daall(ml,ib,'$$DAJUNK$$',inob,invb)
3835	call dapint(ma,ia,ml,ib,jx)
3836	do i=1,ib
3837	call dacop(ml(i),mb(i))
3838	enddo
3839	call dadal(ml,ib)
3840	else
3841	do i=1,ia
3842	call dapek(ma(i),jj,x(i))
3843	call dapok(ma(i),jj,0.0_dp)
3844	enddo
3845	call dapint(ma,ia,mb,ib,jx)
3846	do i=1,ia
3847	call dapok(ma(i),jj,x(i))
3848	enddo
3849	endif
3850	return
3851	end subroutine dapin
3852	subroutine dapint(ma,ia,mb,ib,jind)
3853	implicit none
3854	! **********************************
3855	!
3856	! THIS SUBROUTINE PERFORMS A PARTIAL INVERSION OF THE ROWS MARKED WITH
3857	! NONZERO ENTRIES IN JJ OF THE MATRIX A. THE RESULT IS STORED IN B.
3858	!
3859	!-----------------------------------------------------------------------------
3860	!
3861	integer i,ia,ib,illa,ilma,inoa,inva,ipoa,k
3862	integer,dimension(lnv)::jj,mn,mi,me
3863	integer,dimension(:)::ma,mb,jind
3864	!
3865	if((.not.C_%STABLE_DA)) then
3866	if(c_%watch_user) then
3867	write(6,*) "big problem in dabnew ", sqrt(crash)
3868	endif
3869	return
3870	endif
3871	call dainf(ma(1),inoa,inva,ipoa,ilma,illa)
3872	if((.not.C_%STABLE_DA)) then
3873	if(c_%watch_user) then
3874	write(6,*) "big problem in dabnew ", sqrt(crash)
3875	endif
3876	return
3877	endif
3878	!
3879	do i=1,ia
3880	mn(i) = 0
3881	mi(i) = 0
3882	me(i) = 0
3883	enddo
3884	!
3885	call daall(mn,ia,'$$PIN1 $$',inoa,inva)
3886	call daall(mi,ia,'$$PIN2 $$',inoa,inva)
3887	call daall(me,ia,'$$PIN3 $$',inoa,inva)
3888	!
3889	do i=1,ia
3890	do k=1,nvmax
3891	jj(k) = 0
3892	enddo
3893	jj(i) = 1
3894	call dapok(me(i),jj,1.0_dp)
3895	enddo
3896	!
3897	do i=1,ia
3898	call dacop(ma(i),mn(i))
3899	if(jind(i).eq.0) call dacop(me(i),mn(i))
3900	enddo
3901	!
3902	call dainv(mn,ia,mi,ia)
3903	!
3904	do i=1,ia
3905	if(jind(i).eq.0) call dacop(ma(i),me(i))
3906	enddo
3907	!
3908	call dacct(me,ia,mi,ia,mb,ib)
3909	!
3910	call dadal(me,ia)
3911	call dadal(mi,ia)
3912	call dadal(mn,ia)
3913	!
3914	return
3915	end subroutine dapint
3916	!
3917	subroutine dader(idif,ina,inc)
3918	! subroutine dader(idif,ina,inc)
3919	implicit none
3920	! ******************************
3921	!
3922	! THIS SUBROUTINE COMPUTES THE DERIVATIVE WITH RESPECT TO VARIABLE I
3923	! OF THE VECTOR A AND STORES THE RESULT IN C.
3924	!
3925	!-----------------------------------------------------------------------------
3926	!
3927	integer idif,illc,ilmc,ina,inc,incc,inoc,invc,ipoc
3928	!
3929	if((.not.C_%STABLE_DA)) then
3930	if(c_%watch_user) then
3931	write(6,*) "big problem in dabnew ", sqrt(crash)
3932	endif
3933	return
3934	endif
3935	if(ina.eq.inc) then
3936	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
3937	if((.not.C_%STABLE_DA)) then
3938	if(c_%watch_user) then
3939	write(6,*) "big problem in dabnew ", sqrt(crash)
3940	endif
3941	return
3942	endif
3943	incc=0
3944	call daall1(incc,'$$DAJUNK$$',inoc,invc)
3945	call dadert(idif,ina,incc)
3946	call dacop(incc,inc)
3947	call dadal1(incc)
3948	else
3949	call dadert(idif,ina,inc)
3950	endif
3951	return
3952	end subroutine dader
3953	subroutine dadert(idif,ina,inc)
3954	implicit none
3955	! ******************************
3956	!
3957	! THIS SUBROUTINE COMPUTES THE DERIVATIVE WITH RESPECT TO VARIABLE I
3958	! OF THE VECTOR A AND STORES THE RESULT IN C.
3959	!
3960	!-----------------------------------------------------------------------------
3961	!
3962	integer i,ibase,ic,ider1,ider1s,ider2,ider2s,idif,iee,ifac,illa,illc,&
3963	ilma,ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses
3964	integer,dimension(lnv)::jd
3965	real(dp) rr,x,xdivi
3966	!
3967	if((.not.C_%STABLE_DA)) then
3968	if(c_%watch_user) then
3969	write(6,*) "big problem in dabnew ", sqrt(crash)
3970	endif
3971	return
3972	endif
3973	call dainf(ina,inoa,inva,ipoa,ilma,illa)
3974	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
3975	if((.not.C_%STABLE_DA)) then
3976	if(c_%watch_user) then
3977	write(6,*) "big problem in dabnew ", sqrt(crash)
3978	endif
3979	return
3980	endif
3981	!
3982	!
3983	if(nomax.eq.1) then
3984	! PRINT*,'ERROR, DADER CALLED WITH NOMAX = 1'
3985	! call dadeb !(31,'ERR DADER1',1)
3986	! stop 666
3987	do i=1,lnv
3988	jd(i)=0
3989	enddo
3990	jd(idif)=1
3991	call dapek(ina,jd,rr)
3992	call dacon(inc,rr)
3993	return
3994	endif
3995	!
3996	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC)
3997	!
3998	ibase = nomax + 1
3999	!
4000	if(idif.gt.(nvmax+1)/2) then
4001	ider1 = 0
4002	ider1s = 0
4003	ider2 = idif-(nvmax+1)/2
4004	ider2s = 1
4005	do jj=1,ider2-1
4006	ider2s = ider2s*ibase
4007	enddo
4008	xdivi = ider2s*ibase
4009	else
4010	ider1 = idif
4011	ider1s = 1
4012	do jj=1,ider1-1
4013	ider1s = ider1s*ibase
4014	enddo
4015	ider2 = 0
4016	ider2s = 0
4017	xdivi = ider1s*ibase
4018	endif
4019	!
4020	ibase = nomax+1
4021	!
4022	ic = ipoc-1
4023	!
4024	do i=ipoa,ipoa+illa-1
4025	!
4026	if(ider1.eq.0) then
4027	iee = i_2(i)
4028	else
4029	iee = i_1(i)
4030	endif
4031	!
4032	x = iee/xdivi
4033	ifac = int(ibase*(x-int(x+epsmac)+epsmac))
4034	!
4035	if(ifac.eq.0) goto 100
4036	!
4037	ic = ic + 1
4038	cc(ic) = cc(i)*ifac
4039	i_1(ic) = i_1(i) - ider1s
4040	i_2(ic) = i_2(i) - ider2s
4041	!
4042	100 continue
4043	enddo
4044	!
4045	idall(inc) = ic - ipoc + 1
4046	if(idall(inc).gt.idalm(inc)) then
4047	write(line,'(a15)') 'ERROR IN DADER '
4048	ipause=mypauses(35,line)
4049	call dadeb !(31,'ERR DADER2',1)
4050	endif
4051	!
4052	return
4053	end subroutine dadert
4054	!
4055	!
4056	subroutine dacfuR(ina,fun,inc)
4057	implicit none
4058	! *****************************
4059	!
4060	! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
4061	! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
4062	! RESULT IN C
4063	!
4064	!-----------------------------------------------------------------------------
4065	!
4066	integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
4067	complex(dp),external::fun
4068	!
4069	if((.not.C_%STABLE_DA)) then
4070	if(c_%watch_user) then
4071	write(6,*) "big problem in dabnew ", sqrt(crash)
4072	endif
4073	return
4074	endif
4075	if(ina.eq.inc) then
4076	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
4077	if((.not.C_%STABLE_DA)) then
4078	if(c_%watch_user) then
4079	write(6,*) "big problem in dabnew ", sqrt(crash)
4080	endif
4081	return
4082	endif
4083	incc=0
4084	call daall1(incc,'$$DAJUNK$$',inoc,invc)
4085	call dacfuRt(ina,fun,incc)
4086	call dacop(incc,inc)
4087	call dadal1(incc)
4088	else
4089	call dacfuRt(ina,fun,inc)
4090	endif
4091	return
4092	end subroutine dacfuR
4093	subroutine dacfuRt(ina,fun,inc)
4094	implicit none
4095	! external fun
4096	! *****************************
4097	!
4098	! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
4099	! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
4100	! RESULT IN C
4101	!
4102	!-----------------------------------------------------------------------------
4103	!
4104	integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,&
4105	invc,ipoa,ipoc,ipause,mypauses
4106	integer,dimension(lnv)::j
4107	real(dp) cfac,rr
4108	!
4109	interface
4110	! complex(kind(one)) function fun(abc)
4111	function fun(abc)
4112	use precision_constants
4113	implicit none
4114	complex(dp) fun
4115	integer,dimension(:)::abc
4116	end function fun
4117	end interface
4118	!
4119	if((.not.C_%STABLE_DA)) then
4120	if(c_%watch_user) then
4121	write(6,*) "big problem in dabnew ", sqrt(crash)
4122	endif
4123	return
4124	endif
4125	call dainf(ina,inoa,inva,ipoa,ilma,illa)
4126	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
4127	if((.not.C_%STABLE_DA)) then
4128	if(c_%watch_user) then
4129	write(6,*) "big problem in dabnew ", sqrt(crash)
4130	endif
4131	return
4132	endif
4133	!
4134	!
4135	if(nomax.eq.1) then
4136	do i=1,lnv
4137	j(i)=0
4138	enddo
4139	call dapek(ina,j,rr)
4140	cfac = REAL(fun(j),kind=DP)
4141	rr=cfac*rr
4142	call dapok(inc,j,rr)
4143	do i=1,lnv
4144	j(i)=1
4145	call dapek(ina,j,rr)
4146	cfac = REAL(fun(j),kind=DP)
4147	rr=cfac*rr
4148	call dapok(inc,j,rr)
4149	j(i)=0
4150	enddo
4151	! PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1'
4152	! call dadeb !(31,'ERR DACFU ',1)
4153	! stop 667
4154	return
4155	endif
4156	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC)
4157	!
4158	ic = ipoc - 1
4159	!
4160	do ia=ipoa,ipoa+illa-1
4161	!
4162	call dancd(i_1(ia),i_2(ia),j)
4163	cfac = REAL(fun(j),kind=DP)
4164	! IF(abs(CFAC).LT.EPS) GOTO 100
4165	! IF(abs(CFAC*CC(IA)).LT.EPS) GOTO 100
4166	if(abs(cfac*cc(ia)).lt.eps.or.abs(cc(ia)).lt.eps) goto 100
4167	!
4168	ic = ic + 1
4169	cc(ic) = cc(ia)*cfac
4170	i_1(ic) = i_1(ia)
4171	i_2(ic) = i_2(ia)
4172	!
4173	100 continue
4174	enddo
4175	!
4176	idall(inc) = ic - ipoc + 1
4177	if(idall(inc).gt.idalm(inc)) then
4178	write(line,'(a15)') 'ERROR IN DACFU '
4179	ipause=mypauses(36,line)
4180	call dadeb !(31,'ERR DACFU ',1)
4181	endif
4182	!
4183	return
4184	end subroutine dacfuRt
4185	!
4186	subroutine dacfu(ina,fun,inc)
4187	implicit none
4188	! *****************************
4189	!
4190	! THIS SUBROUTINE APPLIES THE EXTERNAL real(dp) FUNCTION
4191	! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
4192	! RESULT IN C
4193	!
4194	!-----------------------------------------------------------------------------
4195	!
4196	integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
4197	real(dp),external::fun
4198	!
4199	if((.not.C_%STABLE_DA)) then
4200	if(c_%watch_user) then
4201	write(6,*) "big problem in dabnew ", sqrt(crash)
4202	endif
4203	return
4204	endif
4205	if(ina.eq.inc) then
4206	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
4207	if((.not.C_%STABLE_DA)) then
4208	if(c_%watch_user) then
4209	write(6,*) "big problem in dabnew ", sqrt(crash)
4210	endif
4211	return
4212	endif
4213	incc=0
4214	call daall1(incc,'$$DAJUNK$$',inoc,invc)
4215	call dacfut(ina,fun,incc)
4216	call dacop(incc,inc)
4217	call dadal1(incc)
4218	else
4219	call dacfut(ina,fun,inc)
4220	endif
4221	return
4222	end subroutine dacfu
4223	subroutine dacfuI(ina,fun,inc)
4224	implicit none
4225	! *****************************
4226	!
4227	! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
4228	! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
4229	! RESULT IN C
4230	!
4231	!-----------------------------------------------------------------------------
4232	!
4233	integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
4234	complex(dp),external::fun
4235	!
4236	if((.not.C_%STABLE_DA)) then
4237	if(c_%watch_user) then
4238	write(6,*) "big problem in dabnew ", sqrt(crash)
4239	endif
4240	return
4241	endif
4242	if(ina.eq.inc) then
4243	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
4244	if((.not.C_%STABLE_DA)) then
4245	if(c_%watch_user) then
4246	write(6,*) "big problem in dabnew ", sqrt(crash)
4247	endif
4248	return
4249	endif
4250	incc=0
4251	call daall1(incc,'$$DAJUNK$$',inoc,invc)
4252	call dacfuIt(ina,fun,incc)
4253	call dacop(incc,inc)
4254	call dadal1(incc)
4255	else
4256	call dacfuIt(ina,fun,inc)
4257	endif
4258	return
4259	end subroutine dacfuI
4260	subroutine dacfuIt(ina,fun,inc)
4261	implicit none
4262	! external fun
4263	! *****************************
4264	!
4265	! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
4266	! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
4267	! RESULT IN C
4268	!
4269	!-----------------------------------------------------------------------------
4270	!
4271	integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,&
4272	invc,ipoa,ipoc,ipause,mypauses
4273	integer,dimension(lnv)::j
4274	real(dp) cfac,rr
4275	!
4276	interface
4277	! complex(kind(one)) function fun(abc)
4278	function fun(abc)
4279	use precision_constants
4280	implicit none
4281	complex(dp) fun
4282	integer,dimension(:)::abc
4283	end function fun
4284	end interface
4285	!
4286	if((.not.C_%STABLE_DA)) then
4287	if(c_%watch_user) then
4288	write(6,*) "big problem in dabnew ", sqrt(crash)
4289	endif
4290	return
4291	endif
4292	call dainf(ina,inoa,inva,ipoa,ilma,illa)
4293	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
4294	if((.not.C_%STABLE_DA)) then
4295	if(c_%watch_user) then
4296	write(6,*) "big problem in dabnew ", sqrt(crash)
4297	endif
4298	return
4299	endif
4300	!
4301	!
4302	if(nomax.eq.1) then
4303	do i=1,lnv
4304	j(i)=0
4305	enddo
4306	call dapek(ina,j,rr)
4307	cfac = AIMAG(fun(j))
4308	rr=cfac*rr
4309	call dapok(inc,j,rr)
4310	do i=1,lnv
4311	j(i)=1
4312	call dapek(ina,j,rr)
4313	cfac = AIMAG(fun(j))
4314	rr=cfac*rr
4315	call dapok(inc,j,rr)
4316	j(i)=0
4317	enddo
4318	! PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1'
4319	! call dadeb !(31,'ERR DACFU ',1)
4320	! stop 667
4321	return
4322	endif
4323	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC)
4324	!
4325	ic = ipoc - 1
4326	!
4327	do ia=ipoa,ipoa+illa-1
4328	!
4329	call dancd(i_1(ia),i_2(ia),j)
4330	cfac = aimag(fun(j))
4331	! IF(abs(CFAC).LT.EPS) GOTO 100
4332	! IF(abs(CFAC*CC(IA)).LT.EPS) GOTO 100
4333	if(abs(cfac*cc(ia)).lt.eps.or.abs(cc(ia)).lt.eps) goto 100
4334	!
4335	ic = ic + 1
4336	cc(ic) = cc(ia)*cfac
4337	i_1(ic) = i_1(ia)
4338	i_2(ic) = i_2(ia)
4339	!
4340	100 continue
4341	enddo
4342	!
4343	idall(inc) = ic - ipoc + 1
4344	if(idall(inc).gt.idalm(inc)) then
4345	write(line,'(a15)') 'ERROR IN DACFU '
4346	ipause=mypauses(37,line)
4347	call dadeb !(31,'ERR DACFU ',1)
4348	endif
4349	!
4350	return
4351	end subroutine dacfuIt
4352	!
4353	subroutine dacfut(ina,fun,inc)
4354	implicit none
4355	! external fun
4356	! *****************************
4357	!
4358	! THIS SUBROUTINE APPLIES THE EXTERNAL real(dp) FUNCTION
4359	! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
4360	! RESULT IN C
4361	!
4362	!-----------------------------------------------------------------------------
4363	!
4364	integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,&
4365	invc,ipoa,ipoc,ipause,mypauses
4366	integer,dimension(lnv)::j
4367	real(dp) cfac,rr
4368	!
4369	interface
4370	! real(kind(one)) function fun(abc)
4371	function fun(abc)
4372	use precision_constants
4373	implicit none
4374	real(dp) fun
4375	integer,dimension(:)::abc
4376	end function fun
4377	end interface
4378	!
4379	if((.not.C_%STABLE_DA)) then
4380	if(c_%watch_user) then
4381	write(6,*) "big problem in dabnew ", sqrt(crash)
4382	endif
4383	return
4384	endif
4385	call dainf(ina,inoa,inva,ipoa,ilma,illa)
4386	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
4387	if((.not.C_%STABLE_DA)) then
4388	if(c_%watch_user) then
4389	write(6,*) "big problem in dabnew ", sqrt(crash)
4390	endif
4391	return
4392	endif
4393	!
4394	!
4395	if(nomax.eq.1) then
4396	do i=1,lnv
4397	j(i)=0
4398	enddo
4399	call dapek(ina,j,rr)
4400	cfac = fun(j)
4401	rr=cfac*rr
4402	call dapok(inc,j,rr)
4403	do i=1,lnv
4404	j(i)=1
4405	call dapek(ina,j,rr)
4406	cfac = fun(j)
4407	rr=cfac*rr
4408	call dapok(inc,j,rr)
4409	j(i)=0
4410	enddo
4411	! PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1'
4412	! call dadeb !(31,'ERR DACFU ',1)
4413	! stop 667
4414	return
4415	endif
4416	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC)
4417	!
4418	ic = ipoc - 1
4419	!
4420	do ia=ipoa,ipoa+illa-1
4421	!
4422	call dancd(i_1(ia),i_2(ia),j)
4423	cfac = fun(j)
4424	! IF(abs(CFAC).LT.EPS) GOTO 100
4425	! IF(abs(CFAC*CC(IA)).LT.EPS) GOTO 100
4426	if(abs(cfac*cc(ia)).lt.eps.or.abs(cc(ia)).lt.eps) goto 100
4427	!
4428	ic = ic + 1
4429	cc(ic) = cc(ia)*cfac
4430	i_1(ic) = i_1(ia)
4431	i_2(ic) = i_2(ia)
4432	!
4433	100 continue
4434	enddo
4435	!
4436	idall(inc) = ic - ipoc + 1
4437	if(idall(inc).gt.idalm(inc)) then
4438	write(line,'(a15)') 'ERROR IN DACFU '
4439	ipause=mypauses(38,line)
4440	call dadeb !(31,'ERR DACFU ',1)
4441	endif
4442	!
4443	return
4444	end subroutine dacfut
4445	! subroutine GET_C_J(ina,I,C,J)
4446	! subroutine GET_C_J(ina,I,C,J)
4447	! implicit none
4448	!
4449	! INTEGER I,ina
4450	! integer, dimension(lnv)::j
4451	! real(dp) C
4452	!
4453	! C=CC(I)
4454	! call dancd(i_1(I),i_2(I),J)!
4455	! END subroutine GET_C_J
4456	! END subroutine GET_C_J
4457	subroutine dapri(ina,iunit)
4458	implicit none
4459	! ***************************
4460	! Frank
4461	! THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
4462	!
4463	!-----------------------------------------------------------------------------
4464	!
4465	integer i,ii,iii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit,k
4466	integer,dimension(lnv)::j
4467	!
4468	if((.not.C_%STABLE_DA)) then
4469	if(c_%watch_user) then
4470	write(6,*) "big problem in dabnew ", sqrt(crash)
4471	endif
4472	return
4473	endif
4474	if(ina.lt.1.or.ina.gt.nda_dab) then
4475	print*,'ERROR IN DAPRI, INA = ',ina
4476	C_%STABLE_DA=.false.
4477	endif
4478	!
4479	inoa = idano(ina)
4480	inva = idanv(ina)
4481	ipoa = idapo(ina)
4482	ilma = idalm(ina)
4483	illa = idall(ina)
4484	!
4485	write(iunit,'(/1X,A,A,I5,A,I5,A,I5/1X,A/)') daname(ina),', NO =',inoa,', NV =',inva,', INA =',ina,&
4486	'*********************************************'
4487	!
4488	iout = 0
4489	ioa = 0
4490	if(inva.eq.0) then
4491	write(iunit,'(A)') ' I VALUE '
4492	do i = ipoa,ipoa+illa-1
4493	write(iunit,'(I6,2X,G20.13)') i-ipoa, cc(i)
4494	enddo
4495	elseif(nomax.eq.1) then
4496	if(illa.ne.0) write(iunit,'(A)') ' I COEFFICIENT ORDER EXPONENTS'
4497	if(illa.eq.0) write(iunit,'(A)') ' ALL COMPONENTS 0.0_dp '
4498	do i=1,illa
4499	do k=1,inva
4500	j(k)=0
4501	enddo
4502	iout=iout+1
4503	if(i.ne.1) then
4504	j(i-1)=1
4505	ioa=1
4506	endif
4507	write(iunit,'(I6,2X,G20.13,I5,4X,18(2i2,1X))') iout,cc(ipoa+i-1),ioa,(j(iii),iii=1,nvmax)
4508	write(iunit,*) cc(ipoa+i-1)
4509	enddo
4510	else
4511	if(illa.ne.0) write(iunit,'(A)') ' I COEFFICIENT ORDER EXPONENTS'
4512	if(illa.eq.0) write(iunit,'(A)') ' ALL COMPONENTS 0.0_dp '
4513	do ioa = 0,inoa
4514	do ii=ipoa,ipoa+illa-1
4515	if(ieo(ia1(i_1(ii))+ia2(i_2(ii))).ne.ioa) goto 100
4516	call dancd(i_1(ii),i_2(ii),j)
4517	!ETIENNE
4518	if(abs(cc(ii)).gt.eps) then
4519	!ETIENNE
4520	iout = iout+1
4521	write(iunit,'(I6,2X,G20.13,I5,4X,18(2i2,1X))') iout,cc(ii),ioa,(j(iii),iii=1,nvmax)
4522	!ETIENNE
4523	write(iunit,*) cc(ii)
4524	endif
4525	!ETIENNE
4526	!
4527	100 continue
4528	enddo
4529	enddo
4530	!
4531	endif
4532	write(iunit,'(A)') ' '
4533	!
4534	return
4535	end subroutine dapri
4536	subroutine dapri77(ina,iunit)
4537	implicit none
4538	! ***************************
4539	! Etienne
4540	! THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
4541	!
4542	!-----------------------------------------------------------------------------
4543	!
4544	integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit
4545	integer,dimension(lnv)::j
4546	character(10) c10,k10
4547	!
4548	if(iunit.eq.0) return
4549	if((.not.C_%STABLE_DA)) then
4550	if(c_%watch_user) then
4551	write(6,*) "big problem in dabnew ", sqrt(crash)
4552	endif
4553	return
4554	endif
4555	if(ina.lt.1.or.ina.gt.nda_dab) then
4556	write(line,'(a22,i8)') 'ERROR IN DAPRI, INA = ',ina
4557	C_%STABLE_DA=.false.
4558	endif
4559	!
4560	inoa = idano(ina)
4561	inva = idanv(ina)
4562	ipoa = idapo(ina)
4563	ilma = idalm(ina)
4564	illa = idall(ina)
4565	!
4566	write(iunit,'(/1X,A10,A6,I5,A6,I5,A7,I5/1X,A/)') daname(ina),', NO =',inoa,', NV =',inva,', INA =',ina,&
4567	'*********************************************'
4568	!
4569	if(illa.ne.0) write(iunit,'(A)') ' I COEFFICIENT ORDER EXPONENTS'
4570	if(illa.eq.0) write(iunit,'(A)') ' ALL COMPONENTS 0.0_dp '
4571	!
4572	c10=' NO ='
4573	k10=' NV ='
4574	write(iunit,'(A10,I6,A10,I6)') c10,inoa,k10,inva
4575	iout = 0
4576	!
4577	! DO 100 IOA = 0,INOA
4578	do ioa = 0,nocut
4579	do ii=ipoa,ipoa+illa-1
4580	if(nomax.ne.1) then
4581	if(ieo(ia1(i_1(ii))+ia2(i_2(ii))).ne.ioa) goto 100
4582	endif
4583	!ETIENNE
4584	if(abs(cc(ii)).gt.eps) then
4585	!ETIENNE
4586	if(nomax.ne.1) then
4587	call dancd(i_1(ii),i_2(ii),j)
4588	iout = iout+1
4589	else
4590	if(ii.eq.ipoa.and.ioa.eq.1) goto 100
4591	if(ii.gt.ipoa.and.ioa.eq.0) goto 100
4592	do i=1,lnv
4593	j(i)=0
4594	enddo
4595	if(ii.ne.ipoa) j(ii-ipoa)=1
4596	iout = iout+1
4597	endif
4598	!
4599	! WRITE(IUNIT,*) IOA,CC(II),(J(I),I=1,INVA)
4600	if(abs(cc(ii)).gt.eps) then
4601	if(eps.gt.1e-37_dp) then
4602	write(iunit,501) ioa,cc(ii),(j(i),i=1,inva)
4603	else
4604	write(iunit,503) ioa,cc(ii),(j(i),i=1,inva)
4605	endif
4606	endif
4607	501 format(' ', i3,1x,g23.16,1x,100(1x,i2))
4608	503 format(' ', i3,1x,g23.16,1x,100(1x,i2))
4609	502 format(' ', i5,1x,g23.16,1x,100(1x,i2))
4610	endif
4611	!ETIENNE
4612	!
4613	100 continue
4614	enddo
4615	enddo
4616	!
4617	do i=1,lnv
4618	j(i)=0
4619	enddo
4620	if(iout.eq.0) iout=1
4621	write(iunit,502) -iout,0.0_dp,(j(i),i=1,inva)
4622	!
4623	return
4624	end subroutine dapri77
4625	subroutine dashift(ina,inc,ishift)
4626	implicit none
4627	! real(dp) c
4628	! Frank
4629	!
4630	!-----------------------------------------------------------------------------
4631	!
4632	integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,inb,ishift,ich,&
4633	ik,inc,k
4634	integer,dimension(lnv)::j,jd
4635	!
4636	if((.not.C_%STABLE_DA)) then
4637	if(c_%watch_user) then
4638	write(6,*) "big problem in dabnew ", sqrt(crash)
4639	endif
4640	return
4641	endif
4642	inb=0
4643	if(ina.lt.1.or.ina.gt.nda_dab) then
4644	write(line,'(a22,i8)') 'ERROR IN dashift, INA = ',ina
4645	C_%STABLE_DA=.false.
4646	endif
4647	!
4648	inoa = idano(ina)
4649	inva = idanv(ina)
4650	ipoa = idapo(ina)
4651	ilma = idalm(ina)
4652	illa = idall(ina)
4653	call daall1(inb,'$$DAJUNK$$',inoa,inva)
4654	iout = 0
4655	!
4656	! DO 100 IOA = 0,INOA
4657	do ioa = 0,nocut
4658	do ii=ipoa,ipoa+illa-1
4659	if(nomax.ne.1) then
4660	if(ieo(ia1(i_1(ii))+ia2(i_2(ii))).ne.ioa) goto 100
4661	endif
4662	!ETIENNE
4663	if(abs(cc(ii)).gt.eps) then
4664	!ETIENNE
4665	if(nomax.ne.1) then
4666	call dancd(i_1(ii),i_2(ii),j)
4667	iout = iout+1
4668	else
4669	if(ii.eq.ipoa.and.ioa.eq.1) goto 100
4670	if(ii.gt.ipoa.and.ioa.eq.0) goto 100
4671	do i=1,lnv
4672	j(i)=0
4673	enddo
4674	if(ii.ne.ipoa) j(ii-ipoa)=1
4675	iout = iout+1
4676	endif
4677	do k=1,ishift ! put 2004 may
4678	if(j(k)>0 ) then
4679	write(6,*) " trouble in dashift "
4680	stop 888
4681	endif
4682	enddo
4683	!
4684	! WRITE(IUNIT,*) IOA,CC(II),(J(I),I=1,INVA)
4685	if(abs(cc(ii)).gt.eps) then
4686	if(eps.gt.1e-37_dp) then
4687	! write(iunit,501) ioa,cc(ii),(j(i),i=1,inva)
4688	ich=1
4689	do ik=1,ishift
4690	if(j(ik).ne.0) ich=0
4691	enddo
4692	if(ich.eq.1) then
4693	do ik=1,lnv
4694	jd(ik)=0
4695	enddo
4696	do ik=ishift+1,lnv
4697	jd(ik-ishift)=j(ik) !%%%%%%Etienne
4698	enddo
4699	endif
4700	call dapok(inb,jd,cc(ii))
4701	else
4702	! write(iunit,503) ioa,cc(ii),(j(i),i=1,inva)
4703	ich=1
4704	do ik=1,ishift
4705	if(j(ik).ne.0) ich=0
4706	enddo
4707	if(ich.eq.1) then
4708	do ik=1,lnv
4709	jd(ik)=0
4710	enddo
4711	do ik=ishift+1,lnv
4712	jd(ik-ishift)=j(ik) !%%%%%%Etienne
4713	enddo
4714	endif
4715	call dapok(inb,jd,cc(ii))
4716	endif
4717	endif
4718	501 format(' ', i3,1x,g23.16,1x,100(1x,i2))
4719	503 format(' ', i3,1x,g23.16,1x,100(1x,i2))
4720	502 format(' ', i5,1x,g23.16,1x,100(1x,i2))
4721	endif
4722	!ETIENNE
4723	!
4724	100 continue
4725	enddo
4726	enddo
4727	!
4728	do i=1,lnv
4729	j(i)=0
4730	enddo
4731	call dacop(inb,inc)
4732	call dadal1(inb)
4733	!
4734	return
4735	end subroutine dashift
4736	subroutine darea(ina,iunit)
4737	implicit none
4738	! Frank
4739	!-----------------------------------------------------------------------------
4740	!
4741	integer i,ic,iche,ii,ii_1,ii_2,iin,illa,ilma,ina,inoa,inva,io,io1,ipoa,iunit,&
4742	iwarin,iwarno,iwarnv,nno
4743	integer,dimension(lnv)::j
4744	real(dp) c
4745	character(10) c10
4746	if((.not.C_%STABLE_DA)) then
4747	if(c_%watch_user) then
4748	write(6,*) "big problem in dabnew ", sqrt(crash)
4749	endif
4750	return
4751	endif
4752	!
4753	if(ina.lt.1.or.ina.gt.nda_dab) then
4754	C_%STABLE_DA=.false.
4755	print*,'ERROR IN DAREA, INA = ',ina
4756	! X = SQRT(-ONE)
4757	! PRINT*,X
4758	endif
4759	!
4760	inoa = idano(ina)
4761	inva = idanv(ina)
4762	ipoa = idapo(ina)
4763	ilma = idalm(ina)
4764	illa = idall(ina)
4765	!
4766	do i=1,lnv
4767	j(i) = 0
4768	enddo
4769	!
4770	call daclr(1)
4771	call daclr(ina) ! etienne 2008
4772	!
4773	ic = 0
4774	!
4775	iwarno = 0
4776	iwarnv = 0
4777	iwarin = 0
4778	!
4779	read(iunit,'(A10)') c10
4780	read(iunit,'(18X,I4)') nno
4781	read(iunit,'(A10)') c10
4782	read(iunit,'(A10)') c10
4783	read(iunit,'(A10)') c10
4784	!
4785	!
4786	iin = 0
4787	!
4788	10 continue
4789	iin = iin + 1
4790	read(iunit,'(I6,2X,G20.13,I5,4X,18(2i2,1X))') ii,c,io,(j(i),i=1,inva)
4791	!
4792	if(ii.eq.0) goto 20
4793	!ETIENNE
4794	read(iunit,*) c
4795	!ETIENNE
4796	if(ii.ne.iin) then
4797	iwarin = 1
4798	endif
4799	io1 = 0
4800	do i=1,inva
4801	io1 = io1 + j(i)
4802	enddo
4803	!
4804	if(io1.ne.io) then
4805	iwarnv = 1
4806	goto 10
4807	endif
4808	if(io.gt.inoa) then
4809	! IF(IWARNO.EQ.0) PRINT*,'WARNING IN DAREA, FILE ',
4810	! * 'CONTAINS HIGHER ORDERS THAN VECTOR '
4811	iwarno = 1
4812	goto 10
4813	endif
4814	!
4815	if(nomax.ne.1) then
4816	ic = ic + 1
4817	call dadcd(j,ii_1,ii_2)
4818	ic = ia1(ii_1) + ia2(ii_2)
4819	cc(ic) = c
4820	goto 10
4821	else
4822	iche=0
4823	do i=1,inva
4824	if(j(i).eq.1) iche=i
4825	enddo
4826	cc(ipoa+iche)=c
4827	goto 10
4828	endif
4829	!
4830	20 continue
4831	!
4832	if(nomax.ne.1) call dapac(ina)
4833	!
4834	return
4835	end subroutine darea
4836	!FF
4837	!
4838	subroutine darea77(ina,iunit)
4839	implicit none
4840	! ***************************
4841	! Etienne
4842	! THIS SUBROUTINE READS THE DA VECTOR INA FROM UNIT IUNIT.
4843	!
4844	!-----------------------------------------------------------------------------
4845	!
4846	integer i,ic,iche,ii,ii_1,ii_2,iin,illa,ilma,ina,inoa,inva,ipoa,iunit,&
4847	k,nojoh,nvjoh
4848	integer,dimension(lnv)::j
4849	real(dp) c
4850	character(10) c10,k10
4851	!
4852	if((.not.C_%STABLE_DA)) then
4853	if(c_%watch_user) then
4854	write(6,*) "big problem in dabnew ", sqrt(crash)
4855	endif
4856	return
4857	endif
4858	if(ina.lt.1.or.ina.gt.nda_dab) then
4859	write(line,'(a22,i8)') 'ERROR IN darea77, INA = ',ina
4860	C_%STABLE_DA=.false.
4861	endif
4862	!
4863	inoa = idano(ina)
4864	inva = idanv(ina)
4865	ipoa = idapo(ina)
4866	ilma = idalm(ina)
4867	illa = idall(ina)
4868	!
4869	do i=1,lnv
4870	j(i) = 0
4871	enddo
4872	!
4873	call daclr(1)
4874	call daclr(ina) ! etienne 2008
4875	!
4876	!
4877	read(iunit,'(A10)') c10
4878	read(iunit,'(A10)') c10
4879	read(iunit,'(A10)') c10
4880	read(iunit,'(A10)') c10
4881	read(iunit,'(A10)') c10
4882	read(iunit,'(A10,I6,A10,I6)') c10,nojoh,k10,nvjoh
4883	!
4884	iin = 0
4885	!
4886	10 continue
4887	iin = iin + 1
4888	read(iunit,*) ii,c,(j(k),k=1,nvjoh)
4889	if(ii.lt.0) goto 20
4890	do i=inva+1,nvjoh
4891	if(j(i).ne.0) goto 10
4892	enddo
4893	iche=0
4894	do i=1,inva
4895	iche=iche+j(i)
4896	enddo
4897	if(iche.gt.nomax) goto 10
4898	if(nomax.ne.1) then
4899	call dadcd(j,ii_1,ii_2)
4900	ic = ia1(ii_1) + ia2(ii_2)
4901	cc(ic) = c
4902	else
4903	iche=0
4904	do i=1,inva
4905	if(j(i).eq.1) iche=i
4906	enddo
4907	cc(ipoa+iche)=c
4908	endif
4909	goto 10
4910	!
4911	20 continue
4912	!
4913	if(nomax.ne.1) call dapac(ina)
4914	!
4915	return
4916	end subroutine darea77
4917	subroutine dadeb !(iunit,c,istop)
4918	implicit none
4919	! *******************************
4920	!
4921	! THIS SUBROUTINE SERVES AS A DEBUGGING TOOL. IT PRINTS ALL
4922	! NONZERO INFORMATION IN THE COMMON BLOCKS AND ALL DA VECTORS.
4923	!
4924	!-----------------------------------------------------------------------------
4925	!
4926	!integer,dimension(0:1)::i8
4927	!
4928	!etienne
4929	! if(check_da) then
4930	C_%STABLE_DA=.false.
4931	return
4932	! endif
4933	!READ(,) I
4934	!I=SQRT(DBLE(I))
4935	! stop
4936	end subroutine dadeb
4937	!
4938	!
4939	!
4940	!
4941	subroutine dainf(inc,inoc,invc,ipoc,ilmc,illc)
4942	implicit none
4943	! **********************************************
4944	!
4945	! THIS SUBROUTINE SEARCHES THE NUMBER OF DA VECTOR C
4946	! AND RETURS THE INFORMATION IN COMMON DA
4947	!
4948	!-----------------------------------------------------------------------------
4949	!
4950	integer illc,ilmc,inc,inoc,invc,ipoc,ipause,mypauses
4951	!
4952	if(inc.ge.1.and.inc.le.nda_dab) then
4953	inoc = idano(inc)
4954	invc = idanv(inc)
4955	ipoc = idapo(inc)
4956	ilmc = idalm(inc)
4957	illc = idall(inc)
4958	return
4959	endif
4960	!
4961	write(line,'(a26,1x,i8,1x,a11)') 'ERROR IN DAINF, DA VECTOR ',inc,' NOT FOUND '
4962	ipause=mypauses(35,line)
4963	call dadeb !(31,'ERR DAINF ',1)
4964	!
4965	return
4966	end subroutine dainf
4967	!
4968	subroutine dapac(inc)
4969	implicit none
4970	! ************************
4971	!
4972	! THIS SUBROUTINE PACKS THE INFORMATION IN THE SCRATCH VECTOR 1
4973	! INTO THE VECTOR INC. IF LF = 1, THE FILTERING (CF DAMUF) IS
4974	! PERFORMED.
4975	! INVERSE IS DAUNP.
4976	!
4977	!-----------------------------------------------------------------------------
4978	!
4979	integer i,ic,inc,ipoc,ipause,mypauses
4980	real(dp) ccc
4981	!
4982	! call dainf(inc,inoc,invc,ipoc,ilmc,illc)
4983	ipoc = idapo(inc)
4984	!
4985	ic = ipoc - 1
4986	!
4987	do i=1,nmmax
4988	ccc = cc(i)
4989	if(abs(ccc).lt.eps) goto 100
4990	ic = ic + 1
4991	cc(ic) = ccc
4992	i_1(ic) = ie1(i)
4993	i_2(ic) = ie2(i)
4994	100 continue
4995	enddo
4996	!
4997	idall(inc) = ic - ipoc + 1
4998	if(idall(inc).gt.idalm(inc)) then
4999	write(line,'(a15)') 'ERROR IN DAPAC '
5000	ipause=mypauses(35,line)
5001	call dadeb !(31,'ERR DAPAC ',1)
5002	endif
5003	!
5004	return
5005	end subroutine dapac
5006	!
5007	!
5008	subroutine dachk(ina,inoa,inva, inb,inob,invb, inc,inoc,invc)
5009	implicit none
5010	! *************************************************************
5011	!
5012	! THIS SUBROUTINE CHECKS IF THE VECTORS A, B AND C
5013	! HAVE COMPATIBLE ATTRIBUTES
5014	!
5015	!-----------------------------------------------------------------------------
5016	!
5017	integer ierr,ina,inb,inc,inoa,inob,inoc,inva,invb,invc,&
5018	invsum,ipause,mypauses
5019	!
5020	if(lsw.eq.1) return
5021	!
5022	ierr = 0
5023	!
5024	! CASE OF A UNARY OPERATION
5025	! *************************
5026	!
5027	if(inob.eq.-1.and.invb.eq.-1) then
5028	invsum = inva + invc
5029	if(invsum.eq.0) then
5030	if(inoa.gt.inoc) ierr = 1
5031	elseif(invsum.eq.1) then
5032	ierr = 1
5033	else
5034	if(inoa.gt.inoc.or.inva.gt.invc) ierr = 1
5035	endif
5036	if(ierr.eq.1) then
5037	write(line,'(a16,i8,a5,i8,a17,4(1x,i8))') 'ERROR IN DACHK, ',ina,' AND ',inc, &
5038	& ' ARE INCOMPATIBLE',inoa,inva,inoc,invc
5039	ipause=mypauses(35,line)
5040	call dadeb !(31,'ERR DACHK1',1)
5041	endif
5042	!
5043	! CASE OF A BINARY OPERATION
5044	! **************************
5045	!
5046	else
5047	invsum = inva + invb + invc
5048	if(invsum.eq.0) then
5049	if(inoa.gt.inoc.or.inob.gt.inoc) ierr = 1
5050	elseif(invsum.eq.1.or.invsum.eq.2) then
5051	ierr = 1
5052	else
5053	if(inoa.gt.inoc.or.inob.gt.inoc.or.inva.gt.invc.or.invb.gt.invc) ierr = 1
5054	endif
5055	if(ierr.eq.1) then
5056	write(line,'(a16,i8,a1,i8,a5,i8,a17)') 'ERROR IN DACHK, ',ina,',',inb &
5057	& ,' AND ',inc,' ARE INCOMPATIBLE'
5058	ipause=mypauses(35,line)
5059	call dadeb !(31,'ERR DACHK2',1)
5060	endif
5061	endif
5062	!
5063	return
5064	end subroutine dachk
5065	!
5066	subroutine damch(iaa,ia)
5067	implicit none
5068	! ************************
5069	!
5070	! THIS SUBROUTINE CHECKS IF THE IA VECTORS IN THE MATRIX IA HAVE
5071	! IDENTICAL ATTRIBUTES.
5072	!
5073	!-----------------------------------------------------------------------------
5074	!
5075	integer i,ia,illa,ilma,ino1,inoi,inv1,invi,ipoa,ipause,mypauses
5076	integer,dimension(:)::iaa
5077	!
5078	if((.not.C_%STABLE_DA)) then
5079	if(c_%watch_user) then
5080	write(6,*) "big problem in dabnew ", "big problem in dabnew ", sqrt(crash)
5081	endif
5082	return
5083	endif
5084	call dainf(iaa(1),ino1,inv1,ipoa,ilma,illa)
5085	if((.not.C_%STABLE_DA)) then
5086	if(c_%watch_user) then
5087	write(6,*) "big problem in dabnew ", "big problem in dabnew ", sqrt(crash)
5088	endif
5089	return
5090	endif
5091	!
5092	do i=2,ia
5093	call dainf(iaa(i),inoi,invi,ipoa,ilma,illa)
5094	if((.not.C_%STABLE_DA)) then
5095	if(c_%watch_user) then
5096	write(6,*) "big problem in dabnew ", sqrt(crash)
5097	endif
5098	return
5099	endif
5100	if(ino1.ne.inoi.or.inv1.ne.invi) then
5101	write(line,'(a24,i8,a5,i8,a18)') 'ERROR IN DAMCH, VECTORS ',iaa(1), &
5102	& ' AND ',iaa(i),' ARE INCOMPATIBLE '
5103	ipause=mypauses(35,line)
5104	stop
5105	endif
5106	enddo
5107	!
5108	return
5109	end subroutine damch
5110	!
5111	subroutine dadcd(jj,ic1,ic2)
5112	implicit none
5113	! ****************************
5114	!
5115	! THIS SUBROUTINE CODES THE EXPONENTS IN JJ INTO THEIR DA CODES i_1,i_2.
5116	!
5117	!-----------------------------------------------------------------------------
5118	!
5119	integer i,ibase,ic1,ic2,isplit
5120	integer,dimension(lnv)::jj
5121	!
5122	ibase = nomax + 1
5123	isplit = (nvmax+1)/2
5124	ic1 = 0
5125	ic2 = 0
5126	!
5127	do i=nvmax,isplit+1,-1
5128	ic2 = ic2*ibase + jj(i)
5129	enddo
5130	!
5131	do i=isplit,1,-1
5132	ic1 = ic1*ibase + jj(i)
5133	enddo
5134	!
5135	return
5136	end subroutine dadcd
5137	!
5138	subroutine dancd(ic1,ic2,jj)
5139	implicit none
5140	! ****************************
5141	!
5142	! THIS SUBROUTINE ENCODES THE EXPONENTS IN JJ FROM THEIR DA CODES i_1,i_2.
5143	!
5144	!-----------------------------------------------------------------------------
5145	!
5146	integer i,ibase,ic,ic1,ic2,isplit
5147	integer,dimension(:)::jj
5148	real(dp) x
5149	!
5150	ibase = nomax + 1
5151	isplit = (nvmax+1)/2
5152	!
5153	ic = ic1
5154	do i=1,isplit
5155	x = REAL(ic,kind=DP)/REAL(ibase,kind=DP)
5156	ic = int(x+epsmac)
5157	jj(i) = nint(ibase*(x-ic))
5158	enddo
5159	!
5160	ic = ic2
5161	do i=isplit+1,nvmax
5162	x = REAL(ic,kind=DP)/REAL(ibase,kind=DP)
5163	ic = int(x+epsmac)
5164	jj(i) = nint(ibase*(x-ic))
5165	enddo
5166	!
5167	do i=nvmax+1,size(jj) ! 2002.12.2
5168	jj(i) = 0
5169	enddo
5170	!
5171	return
5172	end subroutine dancd
5173	!ETIENNE
5174	subroutine datra(idif,ina,inc)
5175	implicit none
5176	! ******************************
5177	!
5178	! THIS SUBROUTINE COMPUTES THE PSEUDO DERIVATIVE WITH RESPECT TO VARIABLE I
5179	! OF THE VECTOR A AND STORES THE RESULT IN C.
5180	!
5181	! dx^n/dx= x^(n-1)
5182	!-----------------------------------------------------------------------------
5183	!
5184	integer i,ibase,ic,ider1,ider1s,ider2,ider2s,idif,iee,ifac,illa,illc,ilma,&
5185	ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses
5186	real(dp) x,xdivi
5187	!
5188	if((.not.C_%STABLE_DA)) then
5189	if(c_%watch_user) then
5190	write(6,*) "big problem in dabnew ", sqrt(crash)
5191	endif
5192	return
5193	endif
5194	call dainf(ina,inoa,inva,ipoa,ilma,illa)
5195	call dainf(inc,inoc,invc,ipoc,ilmc,illc)
5196	if((.not.C_%STABLE_DA)) then
5197	if(c_%watch_user) then
5198	write(6,*) "big problem in dabnew ", sqrt(crash)
5199	endif
5200	return
5201	endif
5202	!
5203	!
5204	! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC)
5205	!
5206	if(nomax.eq.1) then
5207	call dader(idif,ina,inc)
5208	return
5209	endif
5210	ibase = nomax + 1
5211	!
5212	if(idif.gt.(nvmax+1)/2) then
5213	ider1 = 0
5214	ider1s = 0
5215	ider2 = idif-(nvmax+1)/2
5216	ider2s = 1
5217	do jj=1,ider2-1
5218	ider2s = ider2s*ibase
5219	enddo
5220	xdivi = ider2s*ibase
5221	else
5222	ider1 = idif
5223	ider1s = 1
5224	do jj=1,ider1-1
5225	ider1s = ider1s*ibase
5226	enddo
5227	ider2 = 0
5228	ider2s = 0
5229	xdivi = ider1s*ibase
5230	endif
5231	!
5232	ibase = nomax+1
5233	!
5234	ic = ipoc-1
5235	!
5236	do i=ipoa,ipoa+illa-1
5237	!
5238	if(ider1.eq.0) then
5239	iee = i_2(i)
5240	else
5241	iee = i_1(i)
5242	endif
5243	!
5244	x = iee/xdivi
5245	ifac = int(ibase*(x-int(x+epsmac)+epsmac))
5246	!
5247	if(ifac.eq.0) goto 100
5248	!
5249	!etienne IFAC = INT(IBASE*(X-INT(X)+c_1d_8))
5250	!
5251	!etienne IF(IFAC.EQ.0) GOTO 100
5252	!
5253	ic = ic + 1
5254	cc(ic) = cc(i)
5255	i_1(ic) = i_1(i) - ider1s
5256	i_2(ic) = i_2(i) - ider2s
5257	!
5258	100 continue
5259	enddo
5260	!
5261	idall(inc) = ic - ipoc + 1
5262	if(idall(inc).gt.idalm(inc)) then
5263	write(line,'(a16)') 'ERROR IN DADTRA '
5264	ipause=mypauses(35,line)
5265	call dadeb !(111,'ERR DADTRA',1)
5266	endif
5267	!
5268	return
5269	end subroutine datra
5270	subroutine hash(no1,nv1,jj,ic1,ic2)
5271	implicit none
5272	! ****************************
5273	!
5274	! THIS SUBROUTINE CODES THE EXPONENTS IN JJ INTO THEIR DA CODES i_1,i_2.
5275	!
5276	!-----------------------------------------------------------------------------
5277	!
5278	integer i,ibase,ic1,ic2,isplit,no1,nv1
5279	integer,dimension(:)::jj
5280	!
5281	ibase = no1 + 1
5282	isplit = (nv1+1)/2
5283	ic1 = 0
5284	ic2 = 0
5285	!
5286	do i=nv1,isplit+1,-1
5287	ic2 = ic2*ibase + jj(i)
5288	enddo
5289	!
5290	do i=isplit,1,-1
5291	ic1 = ic1*ibase + jj(i)
5292	enddo
5293	!
5294	return
5295	end subroutine hash
5296	!
5297	subroutine dehash(no1,nv1,ic1,ic2,jj)
5298	implicit none
5299	! ****************************
5300	!
5301	! THIS SUBROUTINE ENCODES THE EXPONENTS IN JJ FROM THEIR DA CODES i_1,i_2.
5302	!
5303	!-----------------------------------------------------------------------------
5304	!
5305	integer i,ibase,ic,ic1,ic2,isplit,no1,nv1
5306	integer,dimension(:)::jj
5307	real(dp) x
5308	!
5309	!frs epsmac=c_1d_7
5310	ibase = no1 + 1
5311	isplit = (nv1+1)/2
5312	!
5313	ic = ic1
5314	do i=1,isplit
5315	x = REAL(ic,kind=DP)/REAL(ibase,kind=DP)
5316	ic = int(x+epsmac)
5317	jj(i) = nint(ibase*(x-ic))
5318	enddo
5319	!
5320	ic = ic2
5321	do i=isplit+1,nv1
5322	x = REAL(ic,kind=DP)/REAL(ibase,kind=DP)
5323	ic = int(x+epsmac)
5324	jj(i) = nint(ibase*(x-ic))
5325	enddo
5326	!
5327	return
5328	end subroutine dehash
5329	subroutine daran(ina,cm,xran)
5330	implicit none
5331	! ************************
5332	!
5333	! THIS SUBROUTINE FILLS THE DA VECTOR A WITH RANDOM ENTRIES.
5334	! FOR CM > 0, THE VECTOR IS FILLED WITH REALS,
5335	! FOR CM < 0, THE VECTOR IS FILLED WITH SINGLE DIGIT INTEGERS
5336	! ABS(CM) IS THE FILLING FACTOR
5337	!
5338	!-----------------------------------------------------------------------------
5339	!
5340	integer i,illa,ilma,ina,inoa,inva,ipoa,ipause,mypauses
5341	real(dp) cm,xran
5342	!
5343	if((.not.C_%STABLE_DA)) then
5344	if(c_%watch_user) then
5345	write(6,*) "big problem in dabnew ", sqrt(crash)
5346	endif
5347	return
5348	endif
5349	call dainf(ina,inoa,inva,ipoa,ilma,illa)
5350	if((.not.C_%STABLE_DA)) then
5351	if(c_%watch_user) then
5352	write(6,*) "big problem in dabnew ", sqrt(crash)
5353	endif
5354	return
5355	endif
5356	!
5357	if(inva.eq.0.or.nomax.eq.1) then
5358	do i=ipoa,ipoa+ilma-1
5359	if(cm.gt.0.0_dp) then
5360	cc(i) = bran(xran)
5361	if(cc(i).gt.cm) cc(i) = 0.0_dp
5362	elseif(cm.lt.0.0_dp) then
5363	cc(i) = int(1+10*bran(xran))
5364	if(cc(i).gt.-10.0_dp*cm) cc(i) = 0.0_dp
5365	endif
5366	enddo
5367	idall(ina) = idalm(ina)
5368	return
5369	endif
5370	!
5371	if(inoa.ne.nomax.or.inva.ne.nvmax) then
5372	line='ERROR IN DARAN, ONLY VECTORS WITH NO = NOMAX AND NV = NVMAX ALLOWED'
5373	ipause=mypauses(31,line)
5374	call dadeb !(31,'ERR DARAN1',1)
5375	endif
5376	!
5377	call daclr(1)
5378	!
5379	do i=1,nmmax
5380	if(cm.gt.0.0_dp) then
5381	cc(i) = bran(xran)
5382	if(cc(i).gt.cm) cc(i) = 0.0_dp
5383	elseif(cm.lt.0.0_dp) then
5384	cc(i) = int(1+10*bran(xran))
5385	if(cc(i).gt.-10.0_dp*cm) cc(i) = 0.0_dp
5386	else
5387	line='ERROR IN ROUTINE DARAN'
5388	ipause=mypauses(31,line)
5389	call dadeb !(31,'ERR DARAN2',1)
5390	endif
5391	enddo
5392	!
5393	call dapac(ina)
5394	!
5395	return
5396	end subroutine daran
5397	!
5398	subroutine dacycle(ina,ipresent,value,illa,j)
5399	implicit none
5400	integer ipause, mypauses
5401	! ***************************
5402	!
5403	! THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
5404	!
5405	!-----------------------------------------------------------------------------
5406	!
5407	integer ii,illa,ilma,ina,inoa,inva,iout,ipoa,ipresent
5408	integer,optional,dimension(:)::j
5409	real(dp) value
5410	!
5411	if(ina.lt.1.or.ina.gt.nda_dab) then
5412	write(line,'(a22,i8)') 'ERROR IN dacycle, INA = ',ina
5413	ipause=mypauses(39,line)
5414	write(6,*) ina
5415	stop
5416	endif
5417	!
5418	inoa = idano(ina)
5419	inva = idanv(ina)
5420	ipoa = idapo(ina)
5421	ilma = idalm(ina)
5422	illa = idall(ina)
5423	if(.not.present(j)) return
5424	iout = 0
5425	if(ipresent.gt.illa.or.ipresent<1) then
5426	write(6,*) ipresent,illa
5427	write(6,*) " error in dacycle "
5428	stop 101
5429	endif
5430	ii=ipresent+ipoa-1
5431	call dancd(i_1(ii),i_2(ii),j)
5432	value=cc(ii)
5433	if(nomax==1) then
5434	j=0
5435	if(ipresent/=1) j(ipresent-1)=1
5436	endif
5437	! ipresent=1+ipresent
5438	return
5439	end subroutine dacycle
5440	!!!! new stuff lingyun
5441	subroutine daclean(ina,value)
5442	implicit none
5443	integer ipause, mypauses
5444	! ***************************
5445	!
5446	! THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
5447	!
5448	!-----------------------------------------------------------------------------
5449	!
5450	integer ii,illa,ilma,ina,inoa,inva,iout,ipoa
5451	real(dp) value
5452	!
5453	if(ina.lt.1.or.ina.gt.nda_dab) then
5454	write(line,'(a22,i8)') 'ERROR IN dacycle, INA = ',ina
5455	ipause=mypauses(39,line)
5456	stop
5457	endif
5458	!
5459	inoa = idano(ina)
5460	inva = idanv(ina)
5461	ipoa = idapo(ina)
5462	ilma = idalm(ina)
5463	illa = idall(ina)
5464	iout = 0
5465	do ii=ipoa,ipoa+ilma-1
5466	if(abs(cc(ii))<value) cc(ii)=0.0_dp
5467	enddo
5468	! call dapac(ina)
5469	return
5470	end subroutine daclean
5471	end module dabnew

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: PSPA/madxPSPA/libs/ptc/src/c_dabnew.f90 @ 430

Download in other formats: