Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

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

Context Navigation

source: PSPA/madxPSPA/src/c_dabnew_berz.f90 @ 430

Download in other formats: