source: PSPA/madxPSPA/libs/ptc/src/j_tpsalie.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: 67.5 KB
Line 
1!The Full Polymorphic Package
2!Copyright (C) Etienne Forest
3
4module tpsalie
5  use tpsa
6  implicit none
7  public
8  private  ASSVEC,ASSMAP,ASSPB,asstaylor,explieflo,expliepb,expflot,exppb
9  private  checkmap,checkpb,checkvec,checktaylor,MAPmatrixr,matrixMAPr,TREEMAP
10  private  TAYLORSMAP,DPEKMAP,DPOKMAP,zeroEQUALMAP,IdentityEQUALMAP
11  private DABSMAP,EQUALMAP,EQUALVEC    !,EQUALMAPVEC,EQUALVECMAP
12  private EQUALvecpb,EQUALpbpb,EQUALpbvec,EQUALpbda,EQUALdapb,CUTORDER,CUTORDERPB,CUTORDERVEC
13  private GETORDERVEC,GETORDERMAP,GETORDERPB,concator,pushtree,concat,pushmatrixr,push1polslow
14  private pushmap
15  private trxflow,trxpb,trxtaylor !,DMULMAPsc,MULMAPsc,IMULMAPsc,DMULVECsc,MULVECsc,IMULVECsc
16  !  private DMULpbsc,MULpbsc,IMULpbsc
17  private scDMULMAP,scMULMAP,scIMULMAP !,scDMULVEC,scMULVEC,scIMULVEC,scDMULpb,scMULpb,scIMULpb
18  private ADDMAP   !,VECMAP,MAPVEC,VECPB,PBVEC
19  private SUBMAP,POWMAP,POWMAP_INV,DAREADTAYLORS,DAREADMAP,DAREADVEC,DAREApb,DAREADTAYLOR
20  PRIVATE DAPRINTTAYLORS,DAPRINTMAP
21  private DAPRINTVEC,DAPRINTTAYLOR,DAPRINTpb,allocmap,allocvec,allocpb,alloctree,allocrad,allocrads
22  private KILLmap,KILLvec,KILLpb,KILLtree,killrad,killrads
23  private A_OPT_damap,k_OPT_damap,A_OPT_vecfield,K_OPT_vecfield,A_OPT_pbfield,K_OPT_pbfield
24  private A_OPT_tree,K_OPT_tree
25  !private push1pol,allocTAYLORS,KILLTAYLORS
26  integer,private::NO,ND,ND2,NP,NDPT,NV
27  logical(lp),private::old
28  !frs real(dp),private::eps=c_1e_9
29  integer::nrmax=400
30  private mul_PBf_t,mul_VECf_t,mul_VECf_MAP,mul_PBf_MAP
31
32
33  private A_OPT_gmap,k_OPT_gmap,allocgmap,KILLgmap,EQUALgMAP,IdentityEQUALgMAP,DAPRINTgMAP,concatorg
34  private assgmap,concatg,DPEKgMAP,DPOKgMAP,gPOWMAP,trxgtaylorc,trxgtaylor,gPOWMAPtpsa,GETORDERgMAP,CUTORDERg
35  private matrixtMAPr,EQUALgMAPdamap
36
37  INTERFACE assignment (=)
38     MODULE PROCEDURE EQUALMAP
39     MODULE PROCEDURE EQUALgMAP
40     MODULE PROCEDURE EQUALgMAPdamap
41     MODULE PROCEDURE MAPTAYLORS
42     MODULE PROCEDURE TAYLORSMAP
43     MODULE PROCEDURE IdentityEQUALMAP
44     MODULE PROCEDURE IdentityEQUALgMAP
45     MODULE PROCEDURE zeroEQUALMAP
46     MODULE PROCEDURE MAPmatrixr
47     MODULE PROCEDURE matrixMAPr
48     MODULE PROCEDURE matrixtMAPr   ! Taylor matrix =  damap
49     !     MODULE PROCEDURE DABSMAP
50     !     MODULE PROCEDURE ABSMAP
51     MODULE PROCEDURE DPEKMAP
52     MODULE PROCEDURE DPOKMAP
53     MODULE PROCEDURE DPEKgMAP
54     MODULE PROCEDURE DPOKgMAP
55     MODULE PROCEDURE EQUALVEC
56     MODULE PROCEDURE EQUALpbpb
57     MODULE PROCEDURE EQUALpbda
58     MODULE PROCEDURE EQUALdapb
59     MODULE PROCEDURE EQUALvecpb
60     MODULE PROCEDURE EQUALpbvec
61     MODULE PROCEDURE TREEMAP
62     !radiation
63     MODULE PROCEDURE radEQUAL
64     MODULE PROCEDURE EQUALrad
65  end  INTERFACE
66
67  INTERFACE OPERATOR (*)
68     !  Move just below  please : here it works
69     MODULE PROCEDURE mul_PBf_MAP  !   Lines to be moved
70     MODULE PROCEDURE mul_PBf_t    !   Lines to be moved
71     MODULE PROCEDURE mul_VECf_t   !   Lines to be moved
72     MODULE PROCEDURE mul_VECf_MAP !   Lines to be moved
73     MODULE PROCEDURE pushmap   ! slow lnv
74     MODULE PROCEDURE pushmatrixr
75     MODULE PROCEDURE push1polslow
76     MODULE PROCEDURE pushtree
77
78     ! DA concatenation
79     MODULE PROCEDURE concat
80     MODULE PROCEDURE concatg
81     MODULE PROCEDURE trxflow
82     MODULE PROCEDURE trxpb
83     MODULE PROCEDURE trxtaylor
84     MODULE PROCEDURE trxgtaylor
85
86
87     MODULE PROCEDURE DMULMAPsc
88     MODULE PROCEDURE MULMAPsc
89     MODULE PROCEDURE IMULMAPsc
90     MODULE PROCEDURE scDMULMAP
91     MODULE PROCEDURE scMULMAP
92     MODULE PROCEDURE scIMULMAP
93     !  Move just below  please
94  END INTERFACE
95
96
97
98  INTERFACE OPERATOR (+)
99     MODULE PROCEDURE ADDMAP
100  END INTERFACE
101
102  INTERFACE OPERATOR (-)
103     MODULE PROCEDURE SUBMAP
104  END INTERFACE
105
106
107  INTERFACE OPERATOR (**)
108     MODULE PROCEDURE POWMAP
109     MODULE PROCEDURE gPOWMAP
110     MODULE PROCEDURE POWMAP_INV
111  END INTERFACE
112
113  INTERFACE OPERATOR (.SUB.)
114     MODULE PROCEDURE GETORDERVEC
115     MODULE PROCEDURE GETORDERMAP
116     MODULE PROCEDURE GETORDERgMAP
117     MODULE PROCEDURE GETORDERPB
118  end  INTERFACE
119
120  INTERFACE OPERATOR (.CUT.)
121     MODULE PROCEDURE CUTORDER
122     MODULE PROCEDURE CUTORDERg
123     MODULE PROCEDURE CUTORDERPB
124     MODULE PROCEDURE CUTORDERVEC
125  END INTERFACE
126
127  INTERFACE OPERATOR (.o.)
128     MODULE PROCEDURE concator
129     MODULE PROCEDURE concatorg
130     MODULE PROCEDURE trxtaylorc
131     MODULE PROCEDURE trxgtaylorc
132     !     MODULE PROCEDURE trxpbc
133     !     MODULE PROCEDURE trxflowc
134  end  INTERFACE
135
136  INTERFACE OPERATOR (.oo.)
137     MODULE PROCEDURE gPOWMAPtpsa
138  end  INTERFACE
139
140  ! i/o
141
142  INTERFACE DAinput
143     MODULE PROCEDURE DAREADTAYLORS
144     MODULE PROCEDURE DAREADMAP
145     MODULE PROCEDURE DAREADVEC
146     MODULE PROCEDURE DAREApb
147     MODULE PROCEDURE DAREADTAYLOR
148  END INTERFACE
149
150  INTERFACE read
151     MODULE PROCEDURE DAREADTAYLORS
152     MODULE PROCEDURE DAREADMAP
153     MODULE PROCEDURE DAREADVEC
154     MODULE PROCEDURE DAREApb
155     MODULE PROCEDURE DAREADTAYLOR
156  END INTERFACE
157
158
159
160  INTERFACE DAprint
161     MODULE PROCEDURE DAPRINTTAYLORS
162     MODULE PROCEDURE DAPRINTMAP
163     MODULE PROCEDURE DAPRINTgMAP
164     MODULE PROCEDURE DAPRINTVEC
165     MODULE PROCEDURE DAPRINTTAYLOR
166     MODULE PROCEDURE DAPRINTpb
167  END INTERFACE
168
169  INTERFACE print
170     MODULE PROCEDURE DAPRINTTAYLORS
171     MODULE PROCEDURE DAPRINTMAP
172     MODULE PROCEDURE DAPRINTgMAP
173     MODULE PROCEDURE DAPRINTVEC
174     MODULE PROCEDURE DAPRINTTAYLOR
175     MODULE PROCEDURE DAPRINTpb
176  END INTERFACE
177
178  ! Exponential of Lie Operators
179
180  INTERFACE texp
181     MODULE PROCEDURE explieflo ! flow on maps
182     MODULE PROCEDURE expliepb  ! pb on maps
183     MODULE PROCEDURE expflot    ! flow on taylor
184     MODULE PROCEDURE exppb     ! pb on taylor
185  END INTERFACE
186
187  INTERFACE exp
188     MODULE PROCEDURE explieflo ! flow on maps
189     MODULE PROCEDURE expliepb  ! pb on maps
190     MODULE PROCEDURE expflot    ! flow on taylor
191     MODULE PROCEDURE exppb     ! pb on taylor
192  END INTERFACE
193
194  INTERFACE full_abs
195     MODULE PROCEDURE DABSMAP
196  END INTERFACE
197
198
199  ! Constructors and Destructors
200
201  INTERFACE alloc
202     MODULE PROCEDURE A_OPT_damap
203     MODULE PROCEDURE A_OPT_gmap
204     MODULE PROCEDURE A_OPT_vecfield
205     MODULE PROCEDURE A_OPT_pbfield
206     MODULE PROCEDURE A_OPT_tree
207     MODULE PROCEDURE allocmap
208     MODULE PROCEDURE allocgmap
209     MODULE PROCEDURE allocvec
210     MODULE PROCEDURE allocpb
211     MODULE PROCEDURE alloctree
212     !radiation
213     MODULE PROCEDURE allocrad
214     MODULE PROCEDURE allocrads
215  END INTERFACE
216
217  INTERFACE allocTPSA
218     MODULE PROCEDURE allocmap
219     MODULE PROCEDURE allocvec
220     MODULE PROCEDURE allocpb
221     MODULE PROCEDURE alloctree
222     !radiation
223     MODULE PROCEDURE allocrad
224     MODULE PROCEDURE allocrads
225  END INTERFACE
226
227  INTERFACE KILL
228     MODULE PROCEDURE k_OPT_damap
229     MODULE PROCEDURE k_OPT_gmap
230     MODULE PROCEDURE k_OPT_vecfield
231     MODULE PROCEDURE k_OPT_pbfield
232     MODULE PROCEDURE k_OPT_tree
233     MODULE PROCEDURE KILLmap
234     MODULE PROCEDURE KILLgmap
235     MODULE PROCEDURE KILLvec
236     MODULE PROCEDURE KILLpb
237     MODULE PROCEDURE KILLtree
238     !radiation
239     MODULE PROCEDURE KILLrad
240     MODULE PROCEDURE KILLrads
241  END INTERFACE
242
243  INTERFACE KILLTPSA
244     MODULE PROCEDURE KILLmap
245     MODULE PROCEDURE KILLvec
246     MODULE PROCEDURE KILLpb
247     MODULE PROCEDURE KILLtree
248     !radiation
249     MODULE PROCEDURE KILLrad
250     MODULE PROCEDURE KILLrads
251  END INTERFACE
252
253  ! Management routines
254
255  INTERFACE ASSDAMAP
256     MODULE PROCEDURE ASSVEC
257     MODULE PROCEDURE ASSMAP
258     MODULE PROCEDURE ASSgMAP
259     MODULE PROCEDURE ASSPB
260     MODULE PROCEDURE asstaylor
261  END INTERFACE
262
263  ! Checking routines
264
265  INTERFACE checkdamap
266     MODULE PROCEDURE checkmap
267     MODULE PROCEDURE checkpb
268     MODULE PROCEDURE checkvec
269     MODULE PROCEDURE checktaylor
270  END INTERFACE
271
272
273contains
274  ! new Poisson stuff
275  FUNCTION mul_PBf_t( S1, S2 )   ! Computes  s1 s2
276    implicit none
277    TYPE (taylor) mul_PBf_t
278    TYPE (PBfield), INTENT (IN) :: S1
279    TYPE (taylor)  , INTENT (IN) :: S2
280    integer localmaster
281    IF(.NOT.C_%STABLE_DA) RETURN
282    localmaster=master
283    call ass(mul_PBf_t)
284
285    mul_PBf_t=0.0_dp
286
287    mul_PBf_t=s1%h.pb.s2
288
289    master=localmaster
290  END FUNCTION mul_PBf_t
291
292  FUNCTION mul_VECf_t( S1, S2 )   ! Computes  s1 s2
293    implicit none
294    TYPE (taylor) mul_VECf_t
295    TYPE (vecfield), INTENT (IN) :: S1
296    TYPE (taylor)  , INTENT (IN) :: S2
297    integer localmaster,i
298    IF(.NOT.C_%STABLE_DA) RETURN
299    localmaster=master
300
301    call ass(mul_VECf_t)
302
303    mul_VECf_t=0.0_dp
304
305    do i=1,c_%nd2
306       mul_VECf_t=mul_VECf_t+s1%v(i)*(s2.d.i)
307    enddo
308    master=localmaster
309  END FUNCTION mul_VECf_t
310
311
312  FUNCTION mul_VECf_MAP( S1, S2 )   ! Computes  s1 s2
313    implicit none
314    TYPE (DAMAP) mul_VECf_MAP
315    TYPE (vecfield), INTENT (IN) :: S1
316    TYPE (DAMAP)  , INTENT (IN) :: S2
317    integer localmaster,i
318    IF(.NOT.C_%STABLE_DA) RETURN
319    localmaster=master
320
321
322    call assDAMAP(mul_VECf_MAP)
323    mul_VECf_MAP=0
324
325    DO i=1,c_%nd2
326       mul_VECf_MAP%V(I)=S1*s2%V(I)
327    ENDDO
328
329
330    master=localmaster
331  END FUNCTION mul_VECf_MAP
332
333  FUNCTION mul_PBf_MAP( S1, S2 )   ! Computes  s1 s2
334    implicit none
335    TYPE (DAMAP) mul_PBf_MAP
336    TYPE (PBfield), INTENT (IN) :: S1
337    TYPE (DAMAP)  , INTENT (IN) :: S2
338    integer localmaster,i
339    IF(.NOT.C_%STABLE_DA) RETURN
340    localmaster=master
341
342
343    call assDAMAP(mul_PBf_MAP)
344    mul_PBf_MAP=0
345
346    DO I=1,C_%ND2
347       mul_PBf_MAP%V(I)=S1*S2%V(I)
348    ENDDO
349    master=localmaster
350  END FUNCTION mul_PBf_MAP
351
352  ! end new Poisson stuff
353
354
355  subroutine set_in_tpsalie( NO1,ND1,ND21,NP1,NDPT1,NV1,log)
356    implicit none
357    integer NO1,ND1,ND21,NP1,NDPT1,NV1
358    logical(lp) log
359    old=log
360    NO=NO1
361    ND=ND1
362    ND2=ND21
363    NP=NP1
364    NDPT=NDPT1
365    NV=NV1
366  end  subroutine set_in_tpsalie
367
368
369
370  SUBROUTINE  alloctree(S1)
371    implicit none
372    type (tree),INTENT(INOUT)::S1
373
374    ! if(old) then
375    call etall(s1%branch%i,nd2)
376    !    else
377    !       call NEWetall(s1%branch%j,nd2)
378    !    endif
379  END SUBROUTINE alloctree
380
381 ! SUBROUTINE  damap_clean(S1,value)
382 !   implicit none
383 !   type (damap),INTENT(INOUT)::S1
384 !   real(dp),INTENT(INOUT)::value
385 !   INTEGER I
386
387!    DO I=1,ND2
388!       CALL taylor_clean(S1%V(I),value)
389!    ENDDO
390
391!  END SUBROUTINE damap_clean
392
393  SUBROUTINE  allocmap(S1)
394    implicit none
395    type (damap),INTENT(INOUT)::S1
396    INTEGER I
397
398    DO I=1,ND2
399       CALL ALLOC(S1%V(I))
400    ENDDO
401    ! if(old) then
402    ! call etall(s1%v%i,nd2)
403    !    else
404    !       call NEWetall(s1%v%j,nd2)
405    !    endif
406
407  END SUBROUTINE allocmap
408
409  SUBROUTINE  allocgmap(S1,n)
410    implicit none
411    type (gmap),INTENT(INOUT)::S1
412    integer, optional :: n
413    integer m
414
415    m=nv
416    if(present(n)) m=n
417    s1%n=m
418
419    ! if(old) then
420    call etall(s1%v%i,s1%n)
421    !    else
422    !       call NEWetall(s1%v%j,s1%n)
423    !    endif
424
425  END SUBROUTINE allocgmap
426
427
428  SUBROUTINE  allocvec(S1)
429    implicit none
430    type (vecfield),INTENT(INOUT)::S1
431
432    s1%ifac=0
433    ! if(old) then
434    call etall(s1%v%i,nd2)
435    !    else
436    !       call NEWetall(s1%v%j,nd2)
437    !    endif
438  END SUBROUTINE allocvec
439
440  SUBROUTINE  allocpb(S1)
441    implicit none
442    type (pbfield),INTENT(INOUT)::S1
443    !     call alloc(s1%h)
444    s1%ifac=0
445    ! if(old) then
446    call etall1(s1%h%i)
447    !    else
448    !       call NEWetall(s1%h%j,1)
449    !    endif
450  END SUBROUTINE allocpb
451
452
453  SUBROUTINE  KILLtree(S1)
454    implicit none
455    type (tree),INTENT(INOUT)::S1
456    ! if(old) then
457    call DADAL(s1%branch%i,nd2)
458    !    else
459    !       call newDADAL(s1%branch%j,nd2)
460    !    endif
461  END SUBROUTINE KILLtree
462
463
464  SUBROUTINE  KILLmap(S1)
465    implicit none
466    type (damap),INTENT(INOUT)::S1
467    INTEGER I
468    ! if(old) then
469    DO I=1,ND2
470       CALL KILL(s1%v(I))
471    ENDDO
472    !    call DADAL(s1%v%i,nd2)
473    !    else
474    !       call newDADAL(s1%v%j,nd2)
475    !    endif
476  END SUBROUTINE KILLmap
477
478  SUBROUTINE  KILLgmap(S1)
479    implicit none
480    type (gmap),INTENT(INOUT)::S1
481    ! if(old) then
482    call DADAL(s1%v%i,s1%n)
483    !    else
484    !       call newDADAL(s1%v%j,s1%n)
485    !    endif
486  END SUBROUTINE KILLgmap
487
488  SUBROUTINE  KILLvec(S1)
489    implicit none
490    type (vecfield),INTENT(INOUT)::S1
491    s1%ifac=0
492    ! if(old) then
493    call DADAL(s1%v%i,nd2)
494    !    else
495    !       call newDADAL(s1%v%j,nd2)
496    !    endif
497  END SUBROUTINE KILLvec
498
499  SUBROUTINE  KILLpb(S1)
500    implicit none
501    type (pbfield),INTENT(INOUT)::S1
502    s1%ifac=0
503    ! if(old) then
504    call DADAL1(s1%h%i)
505    !    else
506    !       call newDADAL(s1%h%j,1)
507    !    endif
508  END SUBROUTINE KILLpb
509
510  SUBROUTINE  A_OPT_damap(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
511    implicit none
512    type (damap),INTENT(INout)::S1,S2
513    type (damap),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
514    call alloc(s1)
515    call alloc(s2)
516    if(present(s3)) call alloc(s3)
517    if(present(s4)) call alloc(s4)
518    if(present(s5)) call alloc(s5)
519    if(present(s6)) call alloc(s6)
520    if(present(s7)) call alloc(s7)
521    if(present(s8)) call alloc(s8)
522    if(present(s9)) call alloc(s9)
523    if(present(s10))call alloc(s10)
524  END SUBROUTINE A_opt_damap
525
526  SUBROUTINE  K_OPT_damap(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
527    implicit none
528    type (damap),INTENT(INout)::S1,S2
529    type (damap),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
530    call KILL(s1)
531    call KILL(s2)
532    if(present(s3)) call KILL(s3)
533    if(present(s4)) call KILL(s4)
534    if(present(s5)) call KILL(s5)
535    if(present(s6)) call KILL(s6)
536    if(present(s7)) call KILL(s7)
537    if(present(s8)) call KILL(s8)
538    if(present(s9)) call KILL(s9)
539    if(present(s10))call KILL(s10)
540  END SUBROUTINE K_OPT_damap
541
542  SUBROUTINE  A_OPT_gmap(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10,n)
543    implicit none
544    type (gmap),INTENT(INout)::S1,S2
545    type (gmap),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
546    integer, optional :: n
547    integer m
548
549    m=nv
550    if(present(n)) m=n
551
552    call alloc(s1,n=m)
553    call alloc(s2,n=m)
554    if(present(s3)) call alloc(s3,n=m)
555    if(present(s4)) call alloc(s4,n=m)
556    if(present(s5)) call alloc(s5,n=m)
557    if(present(s6)) call alloc(s6,n=m)
558    if(present(s7)) call alloc(s7,n=m)
559    if(present(s8)) call alloc(s8,n=m)
560    if(present(s9)) call alloc(s9,n=m)
561    if(present(s10))call alloc(s10,n=m)
562  END SUBROUTINE A_OPT_gmap
563
564  SUBROUTINE  K_OPT_gmap(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
565    implicit none
566    type (gmap),INTENT(INout)::S1,S2
567    type (gmap),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
568    call KILL(s1)
569    call KILL(s2)
570    if(present(s3)) call KILL(s3)
571    if(present(s4)) call KILL(s4)
572    if(present(s5)) call KILL(s5)
573    if(present(s6)) call KILL(s6)
574    if(present(s7)) call KILL(s7)
575    if(present(s8)) call KILL(s8)
576    if(present(s9)) call KILL(s9)
577    if(present(s10))call KILL(s10)
578  END SUBROUTINE K_OPT_gmap
579
580  SUBROUTINE  A_OPT_vecfield(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
581    implicit none
582    type (vecfield),INTENT(INout)::S1,S2
583    type (vecfield),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
584    call alloc(s1)
585    call alloc(s2)
586    if(present(s3)) call alloc(s3)
587    if(present(s4)) call alloc(s4)
588    if(present(s5)) call alloc(s5)
589    if(present(s6)) call alloc(s6)
590    if(present(s7)) call alloc(s7)
591    if(present(s8)) call alloc(s8)
592    if(present(s9)) call alloc(s9)
593    if(present(s10))call alloc(s10)
594  END SUBROUTINE A_OPT_vecfield
595
596  SUBROUTINE  K_OPT_vecfield(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
597    implicit none
598    type (vecfield),INTENT(INout)::S1,S2
599    type (vecfield),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
600    call KILL(s1)
601    call KILL(s2)
602    if(present(s3)) call KILL(s3)
603    if(present(s4)) call KILL(s4)
604    if(present(s5)) call KILL(s5)
605    if(present(s6)) call KILL(s6)
606    if(present(s7)) call KILL(s7)
607    if(present(s8)) call KILL(s8)
608    if(present(s9)) call KILL(s9)
609    if(present(s10))call KILL(s10)
610  END SUBROUTINE K_OPT_vecfield
611
612  SUBROUTINE  A_OPT_pbfield(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
613    implicit none
614    type (pbfield),INTENT(INout)::S1,S2
615    type (pbfield),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
616    call alloc(s1)
617    call alloc(s2)
618    if(present(s3)) call alloc(s3)
619    if(present(s4)) call alloc(s4)
620    if(present(s5)) call alloc(s5)
621    if(present(s6)) call alloc(s6)
622    if(present(s7)) call alloc(s7)
623    if(present(s8)) call alloc(s8)
624    if(present(s9)) call alloc(s9)
625    if(present(s10))call alloc(s10)
626  END SUBROUTINE A_OPT_pbfield
627
628  SUBROUTINE  K_OPT_pbfield(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
629    implicit none
630    type (pbfield),INTENT(INout)::S1,S2
631    type (pbfield),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
632    call KILL(s1)
633    call KILL(s2)
634    if(present(s3)) call KILL(s3)
635    if(present(s4)) call KILL(s4)
636    if(present(s5)) call KILL(s5)
637    if(present(s6)) call KILL(s6)
638    if(present(s7)) call KILL(s7)
639    if(present(s8)) call KILL(s8)
640    if(present(s9)) call KILL(s9)
641    if(present(s10))call KILL(s10)
642  END SUBROUTINE K_OPT_pbfield
643
644  SUBROUTINE  A_OPT_tree(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
645    implicit none
646    type (tree),INTENT(INout)::S1,S2
647    type (tree),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
648    call alloc(s1)
649    call alloc(s2)
650    if(present(s3)) call alloc(s3)
651    if(present(s4)) call alloc(s4)
652    if(present(s5)) call alloc(s5)
653    if(present(s6)) call alloc(s6)
654    if(present(s7)) call alloc(s7)
655    if(present(s8)) call alloc(s8)
656    if(present(s9)) call alloc(s9)
657    if(present(s10))call alloc(s10)
658  END SUBROUTINE A_OPT_tree
659
660  SUBROUTINE  K_OPT_tree(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
661    implicit none
662    type (tree),INTENT(INout)::S1,S2
663    type (tree),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
664    call KILL(s1)
665    call KILL(s2)
666    if(present(s3)) call KILL(s3)
667    if(present(s4)) call KILL(s4)
668    if(present(s5)) call KILL(s5)
669    if(present(s6)) call KILL(s6)
670    if(present(s7)) call KILL(s7)
671    if(present(s8)) call KILL(s8)
672    if(present(s9)) call KILL(s9)
673    if(present(s10))call KILL(s10)
674  END SUBROUTINE K_OPT_tree
675
676
677
678  SUBROUTINE  DAREADTAYLORS(S1,MFILE)
679    implicit none
680    INTEGER,INTENT(in)::MFILE
681    type (TAYLOR),INTENT(INOUT)::S1(NDIM2)
682    INTEGER I
683
684    DO I=1,ND2
685       CALL REA(s1(I),MFILE)
686    ENDDO
687
688  END SUBROUTINE DAREADTAYLORS
689
690  SUBROUTINE  DAREADTAYLOR(S1,MFILE)
691    implicit none
692    INTEGER,INTENT(in)::MFILE
693    type (TAYLOR),INTENT(INOUT)::S1
694
695    CALL REA(s1,MFILE)
696
697  END SUBROUTINE DAREADTAYLOR
698
699  SUBROUTINE  DAREADMAP(S1,MFILE)
700    implicit none
701    INTEGER,INTENT(in)::MFILE
702    type (damap),INTENT(INOUT)::S1
703    INTEGER I
704
705    DO I=1,ND2
706       CALL REA(s1%V(I),MFILE)
707    ENDDO
708
709  END SUBROUTINE DAREADMAP
710
711  SUBROUTINE  DAREADVEC(S1,MFILE)
712    implicit none
713    INTEGER,INTENT(in)::MFILE
714    type (VECFIELD),INTENT(inout)::S1
715    INTEGER I
716
717    read(mfile,*) i
718    s1%ifac=i
719    DO I=1,ND2
720       CALL REA(s1%V(I),MFILE)
721    ENDDO
722
723  END SUBROUTINE DAREADVEC
724
725  SUBROUTINE  DAREAPB(S1,MFILE)
726    implicit none
727    INTEGER,INTENT(IN)::MFILE
728    type (PBFIELD),INTENT(inout)::S1
729    integer i
730    read(mfile,*) i
731    s1%ifac = i
732    CALL REA(s1%H,MFILE)
733
734  END SUBROUTINE DAREAPB
735
736  SUBROUTINE  DAPRINTMAP(S1,MFILE,DEPS)
737    implicit none
738    INTEGER,INTENT(IN)::MFILE
739    type (damap),INTENT(IN)::S1
740    REAL(DP),OPTIONAL,INTENT(INOUT)::DEPS
741    INTEGER I
742
743    DO I=1,ND2
744       CALL PRI(s1%V(I),MFILE,DEPS)
745    ENDDO
746  END SUBROUTINE DAPRINTMAP
747
748  SUBROUTINE  DAPRINTgMAP(S1,MFILE,DEPS)
749    implicit none
750    INTEGER,INTENT(IN)::MFILE
751    type (gmap),INTENT(IN)::S1
752    REAL(DP),OPTIONAL,INTENT(INOUT)::DEPS
753    INTEGER I
754
755    DO I=1,s1%n
756       CALL PRI(s1%V(I),MFILE,DEPS)
757    ENDDO
758  END SUBROUTINE DAPRINTgMAP
759
760  SUBROUTINE  DAPRINTTAYLORS(S1,MFILE,DEPS)
761    implicit none
762    INTEGER,INTENT(IN)::MFILE
763    type (TAYLOR),INTENT(IN)::S1(:)
764    REAL(DP),OPTIONAL,INTENT(INOUT)::DEPS
765    INTEGER I
766
767    DO I=1,size(S1)
768       if(s1(i)%i>0) then
769          if(size(S1)>1) write(MFILE,*) "Taylor #",i
770          CALL PRI(s1(i),MFILE,DEPS)
771       endif
772    ENDDO
773  END SUBROUTINE DAPRINTTAYLORS
774
775  SUBROUTINE  DAPRINTVEC(S1,MFILE,DEPS)
776    implicit none
777    INTEGER,INTENT(IN)::MFILE
778    type (VECFIELD),INTENT(IN)::S1
779    REAL(DP),OPTIONAL,INTENT(INOUT)::DEPS
780    INTEGER I
781
782    write(mfile,*) s1%ifac,' Factorization represented'
783    DO I=1,ND2
784       CALL PRI(s1%V(I),MFILE,DEPS)
785    ENDDO
786
787  END SUBROUTINE DAPRINTVEC
788
789
790  SUBROUTINE  DAPRINTPB(S1,MFILE,DEPS)
791    implicit none
792    INTEGER,INTENT(IN)::MFILE
793    type (PBFIELD),INTENT(IN)::S1
794    REAL(DP),OPTIONAL,INTENT(INOUT)::DEPS
795
796    write(mfile,*) s1%ifac,' Factorization represented'
797    CALL PRI(s1%H,MFILE,DEPS)
798
799  END SUBROUTINE DAPRINTPB
800
801  SUBROUTINE  DAPRINTTAYLOR(S1,MFILE,DEPS)
802    implicit none
803    INTEGER,INTENT(IN)::MFILE
804    type (TAYLOR),INTENT(IN)::S1
805    REAL(DP),OPTIONAL,INTENT(INOUT)::DEPS
806
807    CALL PRI(s1,MFILE,DEPS)
808
809  END SUBROUTINE DAPRINTTAYLOR
810
811  function  DABSMAP(S1)
812    implicit none
813    real(dp) DABSMAP
814    type (damap),INTENT(IN)::S1
815    real(dp) R1
816    INTEGER I
817    IF(.NOT.C_%STABLE_DA) RETURN
818
819    DABSMAP=0.0_dp
820    R1=0.0_dp
821    ! if(old) then
822    if(s1%V(1)%i==0) call crap1("DABSMAP 1")  !call etall(s1%V%i,ND2)
823
824    DO i=1,ND2
825       !          R1=s1%V(I) !2002.10.17
826       R1=full_abs(s1%V(I)) !2002.10.17
827       DABSMAP=DABSMAP+R1
828    ENDDO
829    !    else
830    !       if(.NOT. ASSOCIATED(s1%V(1)%j%r))call crap1("DABSMAP 2")  ! call newetall(s1%V%j,ND2)
831    !
832    !       DO i=1,ND2
833    !          !          R1=s1%V(I)
834    !          R1=full_abs(s1%V(I))
835    !          DABSMAP=DABSMAP+R1
836    !       ENDDO
837    !    endif
838
839  END function DABSMAP
840
841
842  SUBROUTINE  DPEKMAP(S2,S1)
843    implicit none
844    real(dp),INTENT(inOUT),dimension(:)::S2
845    type (damap),INTENT(IN)::S1
846    IF(.NOT.C_%STABLE_DA) RETURN
847    call check_snake
848    ! if(old) then
849    CALL DAPEK0(S1%V%I,S2,nd2)
850    !    else
851    !       CALL newDAPEK0(S1%V%J,S2,nd2)
852    !    endif
853  END SUBROUTINE DPEKMAP
854
855  SUBROUTINE  DPEKgMAP(S2,S1)
856    implicit none
857    real(dp),INTENT(inOUT),dimension(:)::S2
858    type (gmap),INTENT(IN)::S1
859    IF(.NOT.C_%STABLE_DA) RETURN
860    call check_snake
861    ! if(old) then
862    CALL DAPEK0(S1%V%I,S2,s1%n)
863    !    else
864    !       CALL newDAPEK0(S1%V%J,S2,s1%n)
865    !    endif
866  END SUBROUTINE DPEKgMAP
867
868  SUBROUTINE  DPOKMAP(S1,S2)
869    implicit none
870    real(dp),INTENT(IN),dimension(:)::S2
871    type (damap),INTENT(inOUT)::S1
872    IF(.NOT.C_%STABLE_DA) RETURN
873    ! if(old) then
874    if(s1%V(1)%i==0) call crap1("DPOKMAP 1") !call allocw_old(s1%V(1))   !call etall(s1%V%i,ND2)
875    CALL DAPOK0(S1%V%I,S2,nd2)
876    !    else
877    !       if(.NOT. ASSOCIATED(s1%V(1)%j%r))  call crap1("DPOKMAP 2") !  !call allocw_old(s1%V(1))  !call newetall(s1%V%j,ND2)
878    !       CALL NEWDAPOK0(S1%V%J,S2,nd2)
879    !    endif
880  END SUBROUTINE DPOKMAP
881
882  SUBROUTINE  DPOKgMAP(S1,S2)
883    implicit none
884    real(dp),INTENT(IN),dimension(:)::S2
885    type (gmap),INTENT(inOUT)::S1
886    IF(.NOT.C_%STABLE_DA) RETURN
887    ! if(old) then
888    if(s1%V(1)%i==0) call crap1("DPOKMAP 1") !call allocw_old(s1%V(1))   !call etall(s1%V%i,ND2)
889    CALL DAPOK0(S1%V%I,S2,s1%n)
890    !    else
891    !       if(.NOT. ASSOCIATED(s1%V(1)%j%r))  call crap1("DPOKMAP 2") !  !call allocw_old(s1%V(1))  !call newetall(s1%V%j,ND2)
892    !       CALL NEWDAPOK0(S1%V%J,S2,s1%n)
893    !    endif
894  END SUBROUTINE DPOKgMAP
895
896  SUBROUTINE  TREEMAP(S1,S2)
897    implicit none
898    type (damap),INTENT(IN)::S2
899    type (TREE),INTENT(inOUT)::S1
900    integer i
901    IF(.NOT.C_%STABLE_DA) RETURN
902    call check_snake
903    do i=1,nd2
904       call maketree(S2%v(i),S1%branch(i))
905    enddo
906
907  END SUBROUTINE TREEMAP
908
909  SUBROUTINE  matrixtMAPr(S2,S1)
910    implicit none
911    type(taylor),INTENT(inOUT)::S2(:,:)            !(ndim2,ndim2)
912    type (damap),INTENT(IN)::S1
913    integer i,j
914    type(taylor) m(ndim2,ndim2)
915    integer, allocatable :: jl(:)
916
917    IF(.NOT.C_%STABLE_DA) RETURN
918    call check_snake
919
920
921    allocate(JL(nd2))
922    do i=1,nd2
923       do j=1,nd2
924          call alloc(m(i,j))
925       enddo
926    enddo
927    ! if(old) then
928    do i=1,nd2
929       do j=1,nd2
930          JL(j)=1
931          m(i,j)=S1%v(i).par.jl
932          JL(j)=0
933       enddo
934    enddo
935    do i=1,nd2
936       do j=1,nd2
937          s2(i,j)=m(i,j)
938       enddo
939    enddo
940
941    do i=1,nd2
942       do j=1,nd2
943          call kill(m(i,j))
944       enddo
945    enddo
946
947    deallocate(jl)
948
949  END SUBROUTINE matrixtMAPr
950
951
952
953  SUBROUTINE  matrixMAPr(S2,S1)
954    implicit none
955    real(dp),INTENT(inOUT)::S2(:,:)            !(ndim2,ndim2)
956    type (damap),INTENT(IN)::S1
957    integer i,j,JL(lnv)
958    real(dp) m(ndim2,ndim2)
959    IF(.NOT.C_%STABLE_DA) RETURN
960    call check_snake
961
962    do i=1,lnv
963       JL(i)=0
964    enddo
965
966    ! if(old) then
967    do i=1,nd2
968       do j=1,nd2
969          JL(j)=1
970          call dapek(S1%v(i)%i,JL,m(i,j))
971          JL(j)=0
972       enddo
973    enddo
974    do i=1,nd2
975       do j=1,nd2
976          s2(i,j)=m(i,j)
977       enddo
978    enddo
979    !    else
980    !       do i=1,nd2
981    !          do j=1,nd2
982    !             jL(j)=1
983    !             call NEWdapek(S1%v(i)%J,jL,m(i,j))
984    !             jL(j)=0
985    !          enddo
986    !       enddo
987    !       do i=1,nd2
988    !          do j=1,nd2
989    !             s2(i,j)=m(i,j)
990    !          enddo
991    !       enddo
992    !    endif
993  END SUBROUTINE matrixMAPr
994
995  SUBROUTINE  MAPmatrixr(S1,S2)
996    implicit none
997    real(dp),INTENT(in)::S2(:,:)       !    (ndim2,ndim2)
998    type (damap),INTENT(inout)::S1
999    integer i,j,JL(lnv)
1000    IF(.NOT.C_%STABLE_DA) RETURN
1001
1002    do i=1,lnv
1003       JL(i)=0
1004    enddo
1005
1006    s1=0
1007    ! if(old) then
1008    do i=1,nd2
1009       do j=1,nd2
1010          JL(j)=1
1011          call dapok(S1%v(i)%i,JL,s2(i,j))
1012          JL(j)=0
1013       enddo
1014    enddo
1015    !    else
1016    !       do i=1,nd2
1017    !          do j=1,nd2
1018    !             jL(j)=1
1019    !             call NEWdapok(S1%v(i)%J,jL,s2(i,j))
1020    !             jL(j)=0
1021    !         enddo
1022    !       enddo
1023    !    endif
1024  END SUBROUTINE MAPmatrixr
1025
1026
1027  SUBROUTINE  MAPTAYLORS(S2,S1)
1028    implicit none
1029    integer i
1030    type (damap),INTENT(inOUT)::S2
1031    type (TAYLOR),INTENT(IN),dimension(:)::S1
1032    IF(.NOT.C_%STABLE_DA) RETURN
1033    ! if(old) then
1034    if(s2%V(1)%i==0) call crap1("MAPTAYLORS 1")  !call etall(s2%V%i,ND2)
1035    ! CALL DACOPD(S1%I,S2%v%I)
1036    !    else
1037    !       if(.NOT. ASSOCIATED(s2%V(1)%j%r)) call crap1("MAPTAYLORS 2")  !call newetall(s2%V%j,ND2)
1038    !    endif
1039
1040    do i=1,nd2
1041       s2%v(i)=s1(i)
1042    enddo
1043  END SUBROUTINE MAPTAYLORS
1044
1045
1046  SUBROUTINE  TAYLORSMAP(S1,S2)
1047    implicit none
1048    integer i
1049    type (damap),INTENT(IN)::S2
1050    type (TAYLOR),INTENT(inOUT),dimension(:)::S1
1051    IF(.NOT.C_%STABLE_DA) RETURN
1052    call check_snake
1053
1054
1055    ! if(old) then
1056    if(s1(1)%i==0) call crap1("TAYLORSMAP 1")  !call allocw_old(s1(1))
1057    ! CALL DACOPD(S2%v%I,S1%I)
1058    !    else
1059    !       if(.NOT. ASSOCIATED(s1(1)%j%r)) call crap1("TAYLORSMAP 2")  !call allocw_old(s1(1))
1060    !    endif
1061
1062    do i=1,nd2
1063       s1(i)=s2%v(i)
1064    enddo
1065
1066
1067  END SUBROUTINE TAYLORSMAP
1068
1069  SUBROUTINE  EQUALMAP(S2,S1)
1070    implicit none
1071    type (damap),INTENT(inOUT)::S2
1072    type (damap),INTENT(IN)::S1
1073    integer i
1074    IF(.NOT.C_%STABLE_DA) RETURN
1075    call check_snake
1076
1077    ! if(old) then
1078    if(s2%V(1)%i==0) call crap1("EQUALMAP 1")  !call etall(s2%V%i,ND2)
1079    ! CALL DACOPD(S1%V%I,S2%v%I)
1080    !    else
1081    !       if(.NOT. ASSOCIATED(s2%V(1)%J%r)) call crap1("EQUALMAP 2")  !call newetall(s2%V%j,ND2)
1082    !    endif
1083    do i=1,nd2
1084       s2%v(i)=s1%v(i)
1085    enddo
1086  END SUBROUTINE EQUALMAP
1087
1088  SUBROUTINE  EQUALgMAP(S2,S1)
1089    implicit none
1090    type (gmap),INTENT(inOUT)::S2
1091    type (gmap),INTENT(IN)::S1
1092    integer i
1093    IF(.NOT.C_%STABLE_DA) RETURN
1094    call check_snake
1095
1096    do i=1,s1%n
1097       s2%v(i)=s1%v(i)
1098    enddo
1099  END SUBROUTINE EQUALgMAP
1100
1101  SUBROUTINE  EQUALgMAPdamap(S2,S1)
1102    implicit none
1103    type (gmap),INTENT(inOUT)::S2
1104    type (damap),INTENT(IN)::S1
1105    integer i
1106    IF(.NOT.C_%STABLE_DA) RETURN
1107    call check_snake
1108
1109    do i=1,c_%nd2
1110       s2%v(i)=s1%v(i)
1111    enddo
1112  END SUBROUTINE EQUALgMAPdamap
1113
1114
1115
1116  SUBROUTINE  IdentityEQUALMAP(S2,S1)
1117    implicit none
1118    type (damap),INTENT(inOUT)::S2
1119    integer,INTENT(IN)::S1
1120    IF(.NOT.C_%STABLE_DA) RETURN
1121
1122    ! if(old) then
1123    if(s2%V(1)%i==0) call crap1("IdentityEQUALMAP 1") !call etall(s2%V%i,ND2)
1124    IF(S1.EQ.1) CALL etini(S2%V%I)
1125    IF(S1.EQ.0) CALL DACLRD(S2%V%I)
1126    !    else
1127    !       if(.NOT. ASSOCIATED(s2%V(1)%J%r))call crap1("IdentityEQUALMAP 2") ! call newetall(s2%V%j,ND2)
1128    !
1129    !       IF(S1.EQ.1) CALL NEWetini(S2%V%J)
1130    !       IF(S1.EQ.0) CALL NEWDACLRD(S2%V%J)
1131    !    endif
1132  END SUBROUTINE IdentityEQUALMAP
1133
1134  SUBROUTINE  IdentityEQUALgMAP(S2,S1)
1135    implicit none
1136    type (gmap),INTENT(inOUT)::S2
1137    integer,INTENT(IN)::S1
1138    integer i
1139    IF(.NOT.C_%STABLE_DA) RETURN
1140
1141    do i=1,s2%n
1142       if(s1==1) then
1143          s2%v(i)=1.0_dp.mono.i
1144       else
1145          s2%v(i)=0.0_dp
1146
1147       endif
1148    enddo
1149
1150  END SUBROUTINE IdentityEQUALgMAP
1151
1152
1153  SUBROUTINE  zeroEQUALMAP(S2,S1)
1154    implicit none
1155    type (damap),INTENT(inOUT)::S2
1156    real(dp),INTENT(IN)::S1
1157    real(dp) zero_(ndim2)
1158    integer i
1159    IF(.NOT.C_%STABLE_DA) RETURN
1160    do i=1,ndim2
1161       zero_(i)=0.0_dp
1162    enddo
1163    ! if(old) then
1164    if(s2%V(1)%i==0) call crap1("zeroEQUALMAP 1") !call etall(s2%V%i,ND2)
1165    !    else
1166    !       if(.NOT. ASSOCIATED(s2%V(1)%J%r)) call crap1("zeroEQUALMAP 2") !call newetall(s2%V%j,ND2)
1167    !    endif
1168    IF(S1.EQ.0.0_dp) s2=zero_
1169  END SUBROUTINE zeroEQUALMAP
1170
1171
1172  SUBROUTINE  EQUALVEC(S2,S1)
1173    implicit none
1174    type (vecfield),INTENT(inOUT)::S2
1175    type (vecfield),INTENT(IN)::S1
1176    IF(.NOT.C_%STABLE_DA) RETURN
1177    call check_snake
1178    ! if(old) then
1179    if(s2%V(1)%i==0) call crap1("EQUALVEC 1") !call etall(s2%V%i,ND2)
1180    CALL DACOPD(S1%V%I,S2%v%I)
1181    !    else
1182    !       if(.NOT. ASSOCIATED(s2%V(1)%J%r)) call crap1("EQUALVEC 2") !call newetall(s2%V%j,ND2)
1183    !       CALL NEWDACOPD(S1%V%J,S2%v%J)
1184    !    endif
1185    s2%ifac=s1%ifac
1186
1187  END SUBROUTINE EQUALVEC
1188
1189  SUBROUTINE  EQUALvecpb(S2,S1)
1190    implicit none
1191    type (vecfield),INTENT(inOUT)::S2
1192    type (pbfield),INTENT(IN)::S1
1193    IF(.NOT.C_%STABLE_DA) RETURN
1194    call check_snake
1195
1196    ! if(old) then
1197    if(s2%V(1)%i==0) call crap1("EQUALvecpb 1")  !call etall(s2%V%i,ND2)
1198    CALL DIFD(S1%h%i,S2%v%i,-1.0_dp)
1199    !    else
1200    !       if(.NOT. ASSOCIATED(s2%V(1)%J%r)) call crap1("EQUALvecpb 2")  !call newetall(s2%V%j,ND2)
1201    !       CALL NEWDIFD(S1%h%J,S2%v%J,-one)
1202    !    endif
1203    s2%ifac=s1%ifac
1204
1205  END SUBROUTINE EQUALvecpb
1206
1207  SUBROUTINE  EQUALpbvec(S2,S1)
1208    implicit none
1209    type (pbfield),INTENT(inOUT)::S2
1210    type (vecfield),INTENT(IN)::S1
1211
1212    IF(.NOT.C_%STABLE_DA) RETURN
1213    call check_snake
1214
1215    ! if(old) then
1216    if(s2%h%i==0) call crap1("EQUALpbvec 1")  !call etall1(s2%h%i)
1217    CALL intd(S1%v%i,s2%h%i,-1.0_dp)
1218    !    else
1219    !       if(.NOT. ASSOCIATED(s2%h%J%r)) call crap1("EQUALpbvec 2")  !call newetall(s2%h%J,1)
1220    !       CALL NEWintd(S1%v%J,s2%h%J,-one)
1221    !    endif
1222    s2%ifac=s1%ifac
1223  END SUBROUTINE EQUALpbvec
1224
1225
1226  SUBROUTINE  EQUALpbpb(S2,S1)
1227    implicit none
1228    type (pbfield),INTENT(inOUT)::S2
1229    type (pbfield),INTENT(IN)::S1
1230    IF(.NOT.C_%STABLE_DA) RETURN
1231    call check_snake
1232    ! if(old) then
1233    if(s2%h%i==0) call crap1("EQUALpbpb 1")  ! call etall1(s2%h%i)
1234    CALL dacop(S1%h%i,S2%h%i)
1235    !    else
1236    !       if(.NOT. ASSOCIATED(s2%h%J%r)) call crap1("EQUALpbpb 2")  !call newetall(s2%h%J,1)
1237    !       CALL NEWdacop(S1%h%J,S2%h%J)
1238    !    endif
1239    s2%ifac=s1%ifac
1240  END SUBROUTINE EQUALpbpb
1241
1242
1243  SUBROUTINE  EQUALpbda(S2,S1)
1244    implicit none
1245    type (pbfield),INTENT(inOUT)::S2
1246    type (taylor),INTENT(IN)::S1
1247    IF(.NOT.C_%STABLE_DA) RETURN
1248    call check_snake
1249
1250    ! if(old) then
1251    if(s2%h%i==0)  call crap1("EQUALpbda 1")  !call allocw_old(s2%h )
1252    !    else
1253    !       if(.NOT. ASSOCIATED(s2%h%J%r))  call crap1("EQUALpbda 2")  !call allocw_old(s2%h )
1254    !    endif
1255    s2%h=s1
1256  END SUBROUTINE EQUALpbda
1257
1258  SUBROUTINE  EQUALdapb(S2,S1)
1259    implicit none
1260    type (taylor),INTENT(inOUT)::S2
1261    type (pbfield),INTENT(IN)::S1
1262
1263    IF(.NOT.C_%STABLE_DA) RETURN
1264    call check_snake
1265
1266    ! if(old) then
1267    if(s2%i==0) call crap1("EQUALdapb 1")  ! call allocw_old(s2)
1268    !      CALL dacop(S1%h%i,S2%i)
1269    !    else
1270    !       if(.NOT. ASSOCIATED(s2%J%r)) call crap1("EQUALdapb 2")  ! call allocw_old(s2)
1271    !    endif
1272    s2=s1%h
1273
1274  END SUBROUTINE EQUALdapb
1275
1276
1277  FUNCTION CUTORDER( S1, S2 )
1278    implicit none
1279    TYPE (damap) CUTORDER
1280    TYPE (damap), INTENT (IN) :: S1
1281    INTEGER, INTENT (IN) :: S2
1282    INTEGER I
1283    integer localmaster
1284    IF(.NOT.C_%STABLE_DA) RETURN
1285    localmaster=master
1286
1287
1288    call checkdamap(s1)
1289    call assdamap(CUTORDER)
1290
1291    DO I=1,ND2
1292       CUTORDER%V(I)=(S1%V(I)).cut.S2
1293    ENDDO
1294
1295    master=localmaster
1296
1297  END FUNCTION CUTORDER
1298
1299  FUNCTION CUTORDERg( S1, S2 )
1300    implicit none
1301    TYPE (gmap) CUTORDERg
1302    TYPE (gmap), INTENT (IN) :: S1
1303    INTEGER, INTENT (IN) :: S2
1304    INTEGER I
1305    integer localmaster
1306    IF(.NOT.C_%STABLE_DA) RETURN
1307    localmaster=master
1308
1309
1310    call assdamap(CUTORDERg)
1311
1312    DO I=1,ND2
1313       CUTORDERg%V(I)=(S1%V(I)).cut.S2
1314    ENDDO
1315
1316    master=localmaster
1317
1318  END FUNCTION CUTORDERg
1319
1320  FUNCTION CUTORDERPB( S1, S2 )
1321    implicit none
1322    TYPE (PBFIELD) CUTORDERPB
1323    TYPE (PBFIELD), INTENT (IN) :: S1
1324    INTEGER, INTENT (IN) :: S2
1325    integer localmaster
1326    IF(.NOT.C_%STABLE_DA) RETURN
1327    localmaster=master
1328
1329
1330    call checkdamap(s1)
1331    call assdamap(CUTORDERPB)
1332
1333    CUTORDERPB%H=(S1%H).cut.S2
1334    master=localmaster
1335    CUTORDERPB%ifac=S1%ifac
1336
1337  END FUNCTION CUTORDERPB
1338
1339  FUNCTION CUTORDERVEC( S1, S2 )
1340    implicit none
1341    TYPE (VECFIELD) CUTORDERVEC
1342    TYPE (VECFIELD), INTENT (IN) :: S1
1343    INTEGER, INTENT (IN) :: S2
1344    INTEGER I
1345    integer localmaster
1346    IF(.NOT.C_%STABLE_DA) RETURN
1347    localmaster=master
1348
1349
1350    call checkdamap(s1)
1351    call assdamap(CUTORDERVEC)
1352
1353    DO I=1,ND2
1354       CUTORDERVEC%V(I)=(S1%V(I)).cut.S2
1355    ENDDO
1356    CUTORDERVEC%ifac=S1%ifac
1357    master=localmaster
1358
1359  END FUNCTION CUTORDERVEC
1360
1361  FUNCTION GETORDERMAP( S1, S2 )
1362    implicit none
1363    TYPE (damap) GETORDERMAP
1364    TYPE (damap), INTENT (IN) :: S1
1365    INTEGER, INTENT (IN) :: S2
1366    INTEGER I
1367    integer localmaster
1368    IF(.NOT.C_%STABLE_DA) RETURN
1369    localmaster=master
1370
1371
1372    call checkdamap(s1)
1373    call assdamap(GETORDERMAP)
1374
1375    DO I=1,ND2
1376       GETORDERMAP%V(I)=(S1%V(I)).SUB.S2
1377    ENDDO
1378
1379    master=localmaster
1380
1381  END FUNCTION GETORDERMAP
1382
1383  FUNCTION GETORDERgMAP( S1, S2 )
1384    implicit none
1385    TYPE (gmap) GETORDERgMAP
1386    TYPE (gmap), INTENT (IN) :: S1
1387    INTEGER, INTENT (IN) :: S2
1388    INTEGER I
1389    integer localmaster
1390    IF(.NOT.C_%STABLE_DA) RETURN
1391    localmaster=master
1392
1393
1394    call assdamap(GETORDERgMAP)
1395
1396    DO I=1,ND2
1397       GETORDERgMAP%V(I)=(S1%V(I)).SUB.S2
1398    ENDDO
1399
1400    master=localmaster
1401
1402  END FUNCTION GETORDERgMAP
1403
1404  FUNCTION GETORDERVEC( S1, S2 )
1405    implicit none
1406    TYPE (VECFIELD) GETORDERVEC
1407    TYPE (VECFIELD), INTENT (IN) :: S1
1408    INTEGER, INTENT (IN) :: S2
1409    INTEGER I
1410    integer localmaster
1411    IF(.NOT.C_%STABLE_DA) RETURN
1412    localmaster=master
1413
1414
1415    call checkdamap(s1)
1416    call assdamap(GETORDERVEC)
1417
1418    DO I=1,ND2
1419       GETORDERVEC%V(I)=(S1%V(I)).SUB.S2
1420    ENDDO
1421
1422    master=localmaster
1423
1424  END FUNCTION GETORDERVEC
1425
1426  FUNCTION GETORDERPB( S1, S2 )
1427    implicit none
1428    TYPE (PBFIELD) GETORDERPB
1429    TYPE (PBFIELD), INTENT (IN) :: S1
1430    INTEGER, INTENT (IN) :: S2
1431    integer localmaster
1432    IF(.NOT.C_%STABLE_DA) RETURN
1433    localmaster=master
1434
1435
1436
1437    call checkdamap(s1)
1438    call assdamap(GETORDERPB)
1439
1440
1441    GETORDERPB%H=(S1%H).SUB.S2
1442
1443    master=localmaster
1444
1445  END FUNCTION GETORDERPB
1446
1447
1448  FUNCTION pushtree( S1, S2 )
1449    implicit none
1450    TYPE (tree), INTENT (IN) :: S1
1451    real(dp), intent(in),dimension(:)::s2
1452    real(dp) pushtree(lnv) ,junk(lnv)
1453    integer i
1454
1455    do i=1,nd2
1456       junk(i)=push1pol( S1%branch(i), S2 )
1457    enddo
1458
1459
1460    do i=1,nd2
1461       pushtree(I)=junk(i)
1462    enddo
1463    do i=nd2+1,size(s2)
1464       pushtree(I)=S2(i)
1465    enddo
1466
1467  END FUNCTION pushtree
1468
1469  FUNCTION pushmatrixr( S1, S2 )
1470    implicit none
1471    real(dp), INTENT (IN) :: S1(ndim2,ndim2)
1472    real(dp), intent(in)  :: s2(ndim2)
1473    real(dp) pushmatrixr(ndim2) ,junk(ndim2)
1474    integer i,j
1475
1476    pushmatrixr=0.0_dp
1477
1478    do i=1,nd2
1479       junk(i)=0.0_dp
1480    enddo
1481    do i=1,nd2
1482       do j=1,nd2
1483          junk(i)=s1(i,j)*s2(j)+junk(i)
1484       enddo
1485    enddo
1486
1487
1488    do i=1,nd2
1489       pushmatrixr(I)=junk(i)
1490    enddo
1491
1492  END FUNCTION pushmatrixr
1493
1494  FUNCTION pushmap( S1, S2 )
1495    implicit none
1496    TYPE (damap), INTENT (IN) :: S1
1497    TYPE (tree)  arbre
1498    real(dp), intent(in),dimension(:)::s2
1499    real(dp) pushmap(lnv),junk(lnv)
1500    integer i
1501
1502    call alloc(arbre)
1503    arbre=s1
1504    do i=1,nv
1505       junk(i)=s2(i)
1506    enddo
1507    pushmap=arbre*junk
1508    call kill(arbre)
1509  END FUNCTION pushmap
1510
1511
1512  FUNCTION push1pol( S1, S2 )
1513    implicit none
1514    TYPE (taylor), INTENT (IN) :: S1
1515    real(dp), intent(in),dimension(:)::s2
1516    real(dp) push1pol
1517
1518    ! if(old) then
1519    call ppush1(S1%i,s2,push1pol)
1520    !    else
1521    !       call newppush1(S1%j,s2,push1pol)
1522    !    endif
1523
1524
1525  END FUNCTION push1pol
1526
1527  FUNCTION push1polslow( S1, S2 )
1528    implicit none
1529    TYPE (taylor), INTENT (IN) :: S1
1530    real(dp), intent(in),dimension(:)::s2
1531    real(dp) push1polslow
1532    TYPE (taylor) t
1533
1534    call alloc(t)
1535
1536    call maketree(s1,t)
1537    ! if(old) then
1538    call ppush1(t%i,s2,push1polslow)
1539    !    else
1540    !       call newppush1(t%j,s2,push1polslow)
1541    !    endif
1542
1543    call kill(t)
1544
1545  END FUNCTION push1polslow
1546
1547
1548
1549
1550
1551
1552  FUNCTION concat(S1,S2)
1553    implicit none
1554    TYPE (damap) concat,t1,t2,tempnew
1555    TYPE (damap), INTENT (IN) :: S1, S2
1556    real(dp) v1(ndim2),zero_(ndim2)
1557    integer i
1558    integer localmaster
1559    IF(.NOT.C_%STABLE_DA) RETURN
1560    localmaster=master
1561
1562    call checkdamap(s1)
1563    call checkdamap(s2)
1564    call assdamap(concat)
1565    call alloc(t1);call alloc(t2);call alloc(tempnew);
1566
1567
1568    do i=1,ndim2
1569       v1(i)=0.0_dp
1570       zero_(i)=0.0_dp
1571    enddo
1572    t1=s1;t2=s2;
1573    v1=s1     ! change oct 2004.10
1574    t1=zero_;t2=zero_;
1575    ! if(old) then
1576    call etcct(t1%v%i,t2%v%i,tempnew%v%i)
1577    call dacopd(tempnew%v%i,concat%v%i)
1578    !    else
1579    !       call NEWetcct(t1%v%J,t2%v%J,tempnew%v%j)
1580    !       call NEWdacopd(tempnew%v%j,concat%v%J)
1581    !    endif
1582    concat=v1
1583
1584    call kill(t1);call kill(t2);call kill(tempnew);
1585    master=localmaster
1586
1587  END FUNCTION concat
1588
1589  FUNCTION concatg(S1,S2)
1590    implicit none
1591    TYPE (gmap) concatg,t1,t2,tempnew
1592    TYPE (gmap), INTENT (IN) :: S1, S2
1593    real(dp) v1(lnv),zero_(lnv)
1594    integer i
1595    integer localmaster
1596    IF(.NOT.C_%STABLE_DA) RETURN
1597    localmaster=master
1598
1599    concatg%n=s1%n
1600    call assdamap(concatg)
1601    call alloc(t1,S1%n);call alloc(t2,S1%n);call alloc(tempnew,S1%n);
1602
1603
1604    v1=0.0_dp
1605    zero_=0.0_dp
1606
1607    t1=s1;t2=s2;
1608    v1=s2
1609    t1=zero_;t2=zero_;
1610    ! if(old) then
1611    call getcct(t1%v%i,t2%v%i,tempnew%v%i,s1%n)
1612    do i=1,s1%n
1613       call dacop(tempnew%v(i)%i,concatg%v(i)%i)
1614    enddo
1615    !    else
1616    !       call NEWgetcct(t1%v%J,t2%v%J,tempnew%v%j,s1%n)
1617    !       do i=1,s1%n
1618    !          call newdacop(tempnew%v(i)%j,concatg%v(i)%j)
1619    !       enddo
1620    !    endif
1621    concatg=v1
1622
1623    call kill(t1);call kill(t2);call kill(tempnew);
1624    master=localmaster
1625
1626  END FUNCTION concatg
1627
1628  FUNCTION concator( S1, S2 )
1629    implicit none
1630    TYPE (damap) concator
1631    TYPE (damap), INTENT (IN) :: S1, S2
1632    TYPE (damap) tempnew
1633    integer localmaster
1634    IF(.NOT.C_%STABLE_DA) RETURN
1635    localmaster=master
1636
1637    call alloc(tempnew)
1638    call checkdamap(s1)
1639    call checkdamap(s2)
1640    call assdamap(concator)
1641
1642    ! if(old) then
1643    call etcct(s1%v%i,s2%v%i,tempnew%v%i)
1644    call dacopd(tempnew%v%i,concator%v%i)
1645    !    else
1646    !       call NEWetcct(s1%v%J,s2%v%J,tempnew%v%j)
1647    !       call NEWdacopd(tempnew%v%j,concator%v%J)
1648    !    endif
1649
1650    master=localmaster
1651    call kill(tempnew)
1652
1653  END FUNCTION concator
1654
1655  FUNCTION concatorg( S1, S2 )
1656    implicit none
1657    TYPE (gmap) concatorg
1658    TYPE (gmap), INTENT (IN) :: S1, S2
1659    TYPE (gmap) tempnew
1660    integer localmaster,i
1661    IF(.NOT.C_%STABLE_DA) RETURN
1662    localmaster=master
1663
1664    call alloc(tempnew)
1665    concatorg%n=s1%n
1666    call assdamap(concatorg)
1667
1668    ! if(old) then
1669    call getcct(s1%v%i,s2%v%i,tempnew%v%i,s1%n)
1670    do i=1,s1%n
1671       call dacop(tempnew%v(i)%i,concatorg%v(i)%i)
1672    enddo
1673    !    else
1674    !       call NEWgetcct(s1%v%J,s2%v%J,tempnew%v%j,s1%n)
1675    !       do i=1,s1%n
1676    !          call newdacop(tempnew%v(i)%j,concatorg%v(i)%j)
1677    !       enddo
1678    !    endif
1679
1680    master=localmaster
1681    call kill(tempnew)
1682
1683  END FUNCTION concatorg
1684
1685
1686  FUNCTION trxflow(S2,S1)
1687    implicit none
1688    TYPE (vecfield) trxflow
1689    TYPE (vecfield), INTENT (IN) :: S1
1690    TYPE (damap), INTENT (IN) ::  S2
1691    type(damap) s22,tempnew
1692    real(dp) zero_(ndim2)
1693    integer i
1694    integer localmaster
1695    IF(.NOT.C_%STABLE_DA) RETURN
1696    localmaster=master
1697
1698
1699    call alloc(s22,tempnew)
1700    call checkdamap(s1)
1701    call checkdamap(s2)
1702    call assdamap(trxflow)
1703
1704    do i=1,nd2
1705       zero_(i)=0.0_dp
1706    enddo
1707    s22=s2
1708    s22=zero_
1709    ! if(old) then
1710    call trxflo(s1%v%i,tempnew%v%i,s22%v%i)
1711    call dacopd(tempnew%v%i,trxflow%v%i)
1712    !    else
1713    !       call NEWtrxflo(s1%v%J,tempnew%v%J,s22%v%J)
1714    !       call NEWdacopd(tempnew%v%J,trxflow%v%J)
1715    !    endif
1716
1717    call kill(s22,tempnew)
1718    master=localmaster
1719    trxflow%IFAC=S1%IFAC
1720  END FUNCTION trxflow
1721
1722
1723  FUNCTION trxpb( S2, S1 )
1724    implicit none
1725    TYPE (pbfield) trxpb
1726    TYPE (pbfield), INTENT (IN) :: S1
1727    TYPE (damap), INTENT (IN) ::  S2
1728    !    TYPE (damap) S22
1729    !    real(dp) zero_(ndim2)
1730    !    integer i
1731    integer localmaster
1732    IF(.NOT.C_%STABLE_DA) RETURN
1733    localmaster=master
1734
1735
1736    !    call alloc(s22)
1737    call checkdamap(s1)
1738    call checkdamap(s2)
1739    call assdamap(trxpb)
1740
1741
1742    trxpb%h=s1%h*s2
1743
1744    master=localmaster
1745    trxpb%ifac=S1%ifac
1746
1747  END FUNCTION trxpb
1748
1749
1750  FUNCTION trxtaylor( S1, S2 )
1751    implicit none
1752    TYPE (taylor) trxtaylor
1753    TYPE (taylor), INTENT (IN) :: S1
1754    TYPE (damap), INTENT (IN) ::  S2
1755    TYPE (damap)  S22
1756    real(dp) zero_(ndim2)
1757    integer i
1758    integer localmaster
1759    IF(.NOT.C_%STABLE_DA) RETURN
1760    localmaster=master
1761
1762
1763    do i=1,nd2
1764       zero_(i)=0.0_dp
1765    enddo
1766
1767
1768    call checkdamap(s1)
1769    call checkdamap(s2)
1770    call assdamap(trxtaylor)
1771
1772    call alloc(s22)
1773
1774    s22=s2
1775    s22=zero_
1776    ! if(old) then
1777    call trx(s1%i,temp,s22%v%i)
1778    call dacop(temp,trxtaylor%i)
1779    !    else
1780    !       call NEWtrx(s1%J,tempL,s22%v%J)
1781    !       call NEWdacop(tempL,trxtaylor%J)
1782    !   endif
1783
1784    call kill(s22)
1785    master=localmaster
1786
1787  END FUNCTION trxtaylor
1788
1789  FUNCTION trxgtaylor( S1, S2 )
1790    implicit none
1791    TYPE (taylor) trxgtaylor
1792    TYPE (taylor), INTENT (IN) :: S1
1793    TYPE (gmap), INTENT (IN) ::  S2
1794    TYPE (gmap)  S22
1795    real(dp),allocatable:: zero_(:)
1796    integer i
1797    integer localmaster
1798    IF(.NOT.C_%STABLE_DA) RETURN
1799    localmaster=master
1800
1801
1802
1803    call assdamap(trxgtaylor)
1804
1805    call alloc(s22,s2%n)
1806
1807    allocate(zero_(S22%n))
1808
1809    do i=1,S22%n
1810       zero_(i)=0.0_dp
1811    enddo
1812
1813
1814    s22=s2
1815    s22=zero_
1816    ! if(old) then
1817    call gtrx(s1%i,temp,s22%v%i,s2%n)
1818    call dacop(temp,trxgtaylor%i)
1819    !    else
1820    !       call NEWgtrx(s1%J,tempL,s22%v%J,s2%n)
1821    !       call NEWdacop(tempL,trxgtaylor%J)
1822    !   endif
1823
1824    call kill(s22)
1825    deallocate(zero_)
1826    master=localmaster
1827
1828  END FUNCTION trxgtaylor
1829
1830  FUNCTION trxtaylorc( S1, S2 )
1831    implicit none
1832    TYPE (taylor) trxtaylorc
1833    TYPE (taylor), INTENT (IN) :: S1
1834    TYPE (damap), INTENT (IN) ::  S2
1835    integer localmaster
1836    IF(.NOT.C_%STABLE_DA) RETURN
1837    localmaster=master
1838
1839
1840    call checkdamap(s1)
1841    call checkdamap(s2)
1842    call assdamap(trxtaylorc)
1843
1844    ! if(old) then
1845    call trx(s1%i,temp,s2%v%i)
1846    call dacop(temp,trxtaylorc%i)
1847    !    else
1848    !       call NEWtrx(s1%J,tempL,s2%v%J)
1849    !       call NEWdacop(tempL,trxtaylorc%J)
1850    !   endif
1851
1852    master=localmaster
1853
1854  END FUNCTION trxtaylorc
1855
1856  FUNCTION trxgtaylorc( S1, S2 )
1857    implicit none
1858    TYPE (taylor) trxgtaylorc
1859    TYPE (taylor), INTENT (IN) :: S1
1860    TYPE (gmap), INTENT (IN) ::  S2
1861    integer localmaster
1862    IF(.NOT.C_%STABLE_DA) RETURN
1863    localmaster=master
1864
1865
1866    call assdamap(trxgtaylorc)
1867
1868    ! if(old) then
1869    call gtrx(s1%i,temp,s2%v%i,s2%n)
1870    call dacop(temp,trxgtaylorc%i)
1871    !    else
1872    !       call NEWgtrx(s1%J,tempL,s2%v%J,s2%n)
1873    !       call NEWdacop(tempL,trxgtaylorc%J)
1874    !   endif
1875
1876    master=localmaster
1877
1878  END FUNCTION trxgtaylorc
1879
1880
1881
1882  FUNCTION texpdf( S1, S2,NRMIN,NRMAX,SCA,IFAC )
1883    implicit none
1884    TYPE (damap) texpdf
1885    TYPE (vecfield), INTENT (IN) :: S1
1886    TYPE (damap), INTENT (IN) ::  S2
1887    integer NRMIN,NRMAX,IFAC
1888    TYPE (damap) tempnew
1889    real(dp) sca
1890    integer localmaster
1891    IF(.NOT.C_%STABLE_DA) RETURN
1892    localmaster=master
1893
1894
1895    call checkdamap(s1)
1896    call checkdamap(s2)
1897    call assdamap(texpdf)
1898    call alloc(tempnew)
1899
1900    ! if(old) then
1901    !write(6,*) NRMIN,NRMAX,SCA,IFAC
1902    call  FACFLOD(s1%v%i,s2%v%i,tempnew%v%i,NRMIN,NRMAX,SCA,IFAC )
1903    call dacopd(tempnew%v%i,texpdf%v%i)
1904    !    else
1905    !       call newFACFLOD(s1%v%j,s2%v%j,tempnew%v%j,NRMIN,NRMAX,SCA,IFAC )
1906    !       call NEWdacopd(tempnew%v%j,texpdf%v%J)
1907    !   endif
1908    master=localmaster
1909    call kill(tempnew)
1910
1911  END FUNCTION texpdf
1912
1913  FUNCTION texpdft( S1, S2,NRMIN,NRMAX,SCA,IFAC )
1914    implicit none
1915    TYPE (taylor) texpdft
1916    TYPE (vecfield), INTENT (IN) :: S1
1917    TYPE (taylor), INTENT (IN) ::  S2
1918    integer NRMIN,NRMAX,IFAC
1919    real(dp) sca
1920    integer localmaster
1921    IF(.NOT.C_%STABLE_DA) RETURN
1922    localmaster=master
1923
1924
1925    call checkdamap(s1)
1926    call checkdamap(s2)
1927    call assdamap(texpdft)
1928
1929    ! if(old) then
1930    !write(6,*) NRMIN,NRMAX,SCA,IFAC
1931    call  FACFLO(s1%v%i,s2%i,temp,NRMIN,NRMAX,SCA,IFAC )
1932    call dacop(temp,texpdft%i)
1933    !    else
1934    !       call newFACFLO(s1%v%j,s2%j,tempL,NRMIN,NRMAX,SCA,IFAC )
1935    !       call NEWdacop(tempL,texpdft%J)
1936    !   endif
1937    master=localmaster
1938
1939  END FUNCTION texpdft
1940
1941
1942  FUNCTION explieflo( S1, S2 )
1943    implicit none
1944    TYPE (damap) explieflo
1945    TYPE (vecfield), INTENT (IN) :: S1
1946    TYPE (damap), INTENT (IN) ::  S2
1947    TYPE (damap) tempnew
1948
1949    integer no1
1950    integer localmaster
1951    IF(.NOT.C_%STABLE_DA) RETURN
1952    localmaster=master
1953
1954    no1=no
1955
1956    call checkdamap(s1)
1957    call checkdamap(s2)
1958    call assdamap(explieflo)
1959    call alloc(tempnew)
1960    if(s1%ifac/=0) then
1961       explieflo=texpdf( S1, S2,2,no1,1.0_dp,s1%ifac )
1962       !         write(6,*) "Improper usage: map is factorized "
1963       !     write(6,*) "Sorry, we will implement in later version "
1964       !     write(6,*) "Use a Dragt-Finn or Reverse-Dragt-Finn scratch variable "
1965    else
1966       ! if(old) then
1967       call EXPFLOD(s1%v%i,s2%v%i,tempnew%v%i,eps_tpsalie,nrmax)
1968       call dacopd(tempnew%v%i,explieflo%v%i)
1969       !       else
1970       !          call NEWEXPFLOD(s1%v%J,s2%v%J,tempnew%v%J,eps_tpsalie,nrmax)
1971       !          call NEWdacopd(tempnew%v%J,explieflo%v%J)
1972       !      endif
1973    endif
1974    master=localmaster
1975    call kill(tempnew)
1976
1977  END FUNCTION explieflo
1978
1979  FUNCTION expflot( S1, S2 )
1980    implicit none
1981    TYPE (taylor) expflot
1982    TYPE (vecfield), INTENT (IN) :: S1
1983    TYPE (taylor), INTENT (IN) ::  S2
1984    integer no1
1985    integer localmaster
1986    IF(.NOT.C_%STABLE_DA) RETURN
1987    localmaster=master
1988
1989
1990    no1=no
1991    call checkdamap(s1)
1992    call checkdamap(s2)
1993    call assdamap(expflot)
1994
1995    if(s1%ifac/=0) then
1996       expflot=texpdft( S1, S2,2,no1,1.0_dp,s1%ifac )
1997       !         write(6,*) "Improper usage: map is factorized "
1998       !    write(6,*) "Sorry, we will implement in later version "
1999       !     write(6,*) "Use a Dragt-Finn or Reverse-Dragt-Finn scratch variable "
2000    else
2001       ! if(old) then
2002       call EXPFLO(s1%v%i,s2%i,temp,eps_tpsalie,nrmax)
2003       call dacop(temp,expflot%i)
2004       !      else
2005       !         call NEWEXPFLO(s1%v%J,s2%J,tempL,eps_tpsalie,nrmax)
2006       !         call NEWdacop(tempL,expflot%J)
2007       !      endif
2008    endif
2009    master=localmaster
2010
2011  END FUNCTION expflot
2012
2013  FUNCTION expliepb( S1, S2 )
2014    implicit none
2015    TYPE (damap) expliepb
2016    TYPE (pbfield), INTENT (IN) :: S1
2017    TYPE (damap), INTENT (IN) ::  S2
2018    type (vecfield) T
2019    integer localmaster
2020    IF(.NOT.C_%STABLE_DA) RETURN
2021    localmaster=master
2022
2023
2024    call alloc(T)
2025
2026    call checkdamap(s1)
2027    call checkdamap(s2)
2028    call assdamap(expliepb)
2029    T=S1
2030    expliepb=explieflo(T,S2)
2031
2032    CALL KILL(T)
2033    master=localmaster
2034
2035  END FUNCTION expliepb
2036
2037  FUNCTION exppb( S1, S2 )
2038    implicit none
2039    TYPE (taylor) exppb
2040    TYPE (pbfield), INTENT (IN) :: S1
2041    TYPE (taylor), INTENT (IN) ::  S2
2042    type (vecfield) T
2043    integer localmaster
2044    IF(.NOT.C_%STABLE_DA) RETURN
2045    localmaster=master
2046
2047
2048    call alloc(T)
2049
2050    call checkdamap(s1)
2051    call checkdamap(s2)
2052    call assdamap(exppb)
2053    T=S1
2054    exppb=expfloT(T,S2)
2055
2056    CALL KILL(T)
2057    master=localmaster
2058
2059  END FUNCTION exppb
2060
2061
2062  FUNCTION ADDMAP( S1, S2 )
2063    implicit none
2064    TYPE (damap) ADDMAP
2065    TYPE (damap), INTENT (IN) :: S1, S2
2066    integer localmaster
2067    IF(.NOT.C_%STABLE_DA) RETURN
2068    localmaster=master
2069
2070
2071    call checkdamap(s1)
2072    call checkdamap(s2)
2073    call assdamap(ADDMAP)
2074
2075    ! if(old) then
2076    call DALIND(s1%v%i,1.0_dp,s2%v%i,1.0_dp,ADDMAP%v%i)
2077    !   else
2078    !      call newDALIND(s1%v%j,one,s2%v%j,one,ADDMAP%v%j)
2079    !   endif
2080
2081    master=localmaster
2082
2083  END FUNCTION ADDMAP
2084
2085  FUNCTION SUBMAP( S1, S2 )
2086    implicit none
2087    TYPE (damap) SUBMAP
2088    TYPE (damap), INTENT (IN) :: S1, S2
2089    integer localmaster
2090    IF(.NOT.C_%STABLE_DA) RETURN
2091    localmaster=master
2092
2093
2094    call checkdamap(s1)
2095    call checkdamap(s2)
2096    call assdamap(SUBMAP)
2097
2098    ! if(old) then
2099    call DALIND(s1%v%i,1.0_dp,s2%v%i,-1.0_dp,SUBMAP%v%i)
2100    !   else
2101    !      call newDALIND(s1%v%j,one,s2%v%j,-one,SUBMAP%v%j)
2102    !   endif
2103
2104    master=localmaster
2105
2106  END FUNCTION SUBMAP
2107
2108  FUNCTION DMULMAPsc( S1, Sc )
2109    implicit none
2110    TYPE (damap)  DMULMAPsc
2111    TYPE (damap), INTENT (IN) :: S1
2112    real(dp), INTENT (IN) :: SC
2113    INTEGER I
2114    integer localmaster
2115    TYPE (damap)  tempnew
2116    IF(.NOT.C_%STABLE_DA) RETURN
2117    localmaster=master
2118
2119
2120    call checkdamap(s1)
2121    call assdamap(DMULMAPsc)
2122    call alloc(tempnew)
2123
2124    ! if(old) then
2125    DO I=1,ND2
2126       CALL DACMU(s1%v(I)%i,SC,tempnew%v(I)%i)
2127    ENDDO
2128    call dacopd(tempnew%v%i,DMULMAPsc%v%i)
2129    !    else
2130    !       DO I=1,ND2
2131    !          CALL NEWDACMU(s1%v(I)%J,SC,tempnew%v(I)%J)
2132    !       ENDDO
2133    !       call NEWdacopd(tempnew%v%J,DMULMAPsc%v%J)
2134    !    endif
2135
2136    master=localmaster
2137    call kill(tempnew)
2138
2139  END FUNCTION DMULMAPsc
2140
2141  FUNCTION MULMAPsc( S1, Sc )
2142    implicit none
2143    TYPE (damap)  MULMAPsc
2144    TYPE (damap), INTENT (IN) :: S1
2145    REAL(SP), INTENT (IN) :: SC
2146    INTEGER I
2147    TYPE (damap)  tempnew
2148    integer localmaster
2149    IF(.NOT.C_%STABLE_DA) RETURN
2150    localmaster=master
2151
2152
2153    call checkdamap(s1)
2154    call assdamap(MULMAPsc)
2155    call alloc(tempnew)
2156
2157    ! if(old) then
2158    DO I=1,ND2
2159       CALL DACMU(s1%v(I)%i,REAL(SC,kind=DP),tempnew%v(I)%i)
2160    ENDDO
2161    call dacopd(tempnew%v%i,MULMAPsc%v%i)
2162    !    else
2163    !       DO I=1,ND2
2164    !          CALL NEWDACMU(s1%v(I)%J,REAL(SC,kind=DP),tempnew%v(I)%j)
2165    !       ENDDO
2166    !       call NEWdacopd(tempnew%v%j,MULMAPsc%v%J)
2167    !    endif
2168    master=localmaster
2169    call kill(tempnew)
2170
2171  END FUNCTION MULMAPsc
2172
2173  FUNCTION IMULMAPsc( S1, Sc )
2174    implicit none
2175    TYPE (damap)  IMULMAPsc
2176    TYPE (damap), INTENT (IN) :: S1
2177    INTEGER, INTENT (IN) :: SC
2178    INTEGER I
2179    TYPE (damap)  tempnew
2180    integer localmaster
2181    IF(.NOT.C_%STABLE_DA) RETURN
2182    localmaster=master
2183
2184
2185    call checkdamap(s1)
2186    call assdamap(IMULMAPsc)
2187    call alloc(tempnew)
2188    ! if(old) then
2189    DO I=1,ND2
2190       CALL DACMU(s1%v(I)%i,REAL(SC,kind=DP),tempnew%v(I)%i)
2191    ENDDO
2192    call dacopd(tempnew%v%i,IMULMAPsc%v%i)
2193    !    else
2194    !       DO I=1,ND2
2195    !          CALL NEWDACMU(s1%v(I)%J,REAL(SC,kind=DP),tempnew%v(I)%j)
2196    !       ENDDO
2197    !       call NEWdacopd(tempnew%v%j,IMULMAPsc%v%J)
2198    !    endif
2199    master=localmaster
2200    call kill(tempnew)
2201
2202  END FUNCTION IMULMAPsc
2203
2204  FUNCTION scDMULMAP( Sc , S1 )
2205    implicit none
2206    TYPE (damap)  scDMULMAP
2207    TYPE (damap), INTENT (IN) :: S1
2208    real(dp), INTENT (IN) :: SC
2209    INTEGER I
2210    TYPE (damap)  tempnew
2211    integer localmaster
2212    IF(.NOT.C_%STABLE_DA) RETURN
2213    localmaster=master
2214
2215    call checkdamap(s1)
2216    call assdamap(scDMULMAP)
2217
2218    call alloc(tempnew)
2219    ! if(old) then
2220    DO I=1,ND2
2221       CALL DACMU(s1%v(I)%i,SC,tempnew%v(I)%i)
2222    ENDDO
2223    call dacopd(tempnew%v%i,scDMULMAP%v%i)
2224    !    else
2225    !       DO I=1,ND2
2226    !          CALL NEWDACMU(s1%v(I)%J,SC,tempnew%v(i)%j)
2227    !       ENDDO
2228    !       call NEWdacopd(tempnew%v%j,scDMULMAP%v%J)
2229    !    endif
2230    master=localmaster
2231    call kill(tempnew)
2232
2233  END FUNCTION scDMULMAP
2234
2235  FUNCTION scMULMAP( Sc , S1 )
2236    implicit none
2237    TYPE (damap)  scMULMAP
2238    TYPE (damap), INTENT (IN) :: S1
2239    REAL(SP), INTENT (IN) :: SC
2240    INTEGER I
2241    integer localmaster
2242    TYPE (damap)  tempnew
2243    IF(.NOT.C_%STABLE_DA) RETURN
2244    localmaster=master
2245
2246
2247    call checkdamap(s1)
2248    call assdamap(scMULMAP)
2249    call alloc(tempnew)
2250
2251    ! if(old) then
2252    DO I=1,ND2
2253       CALL DACMU(s1%v(I)%i,REAL(SC,kind=DP),tempnew%v(I)%i)
2254    ENDDO
2255    call dacopd(tempnew%v%i,scMULMAP%v%i)
2256    !    else
2257    !       DO I=1,ND2
2258    !          CALL NEWDACMU(s1%v(I)%J,REAL(SC,kind=DP),tempnew%v(I)%j)
2259    !       ENDDO
2260    !       call NEWdacopd(tempnew%v%j,scMULMAP%v%J)
2261    !    endif
2262    master=localmaster
2263    call kill(tempnew)
2264
2265  END FUNCTION scMULMAP
2266
2267  FUNCTION scIMULMAP( Sc , S1 )
2268    implicit none
2269    TYPE (damap)  scIMULMAP
2270    TYPE (damap), INTENT (IN) :: S1
2271    INTEGER, INTENT (IN) :: SC
2272    INTEGER I
2273    integer localmaster
2274    TYPE (damap)  tempnew
2275    IF(.NOT.C_%STABLE_DA) RETURN
2276    localmaster=master
2277
2278
2279    call checkdamap(s1)
2280    call assdamap(scIMULMAP)
2281    call alloc(tempnew)
2282
2283    ! if(old) then
2284    DO I=1,ND2
2285       CALL DACMU(s1%v(I)%i,REAL(SC,kind=DP),tempnew%v(I)%i)
2286    ENDDO
2287    call dacopd(tempnew%v%i,scIMULMAP%v%i)
2288    !    else
2289    !       DO I=1,ND2
2290    !          CALL NEWDACMU(s1%v(I)%J,REAL(SC,kind=DP),tempnew%v(I)%j)
2291    !       ENDDO
2292    !       call NEWdacopd(tempnew%v%j,scIMULMAP%v%J)
2293    !    endif
2294    master=localmaster
2295    call kill(tempnew)
2296
2297  END FUNCTION scIMULMAP
2298
2299
2300
2301  FUNCTION POWMAP( S1, R2 )
2302    implicit none
2303    TYPE (damap) POWMAP
2304    TYPE (damap), INTENT (IN) :: S1
2305    INTEGER, INTENT (IN) :: R2
2306    TYPE (damap) S11
2307    INTEGER I,R22
2308    integer localmaster
2309    IF(.NOT.C_%STABLE_DA) RETURN
2310    localmaster=master
2311
2312
2313    call checkdamap(s1)
2314
2315    call assdamap(POWMAP)
2316
2317    call alloc(s11)
2318
2319    s11=1
2320
2321
2322    R22=IABS(R2)
2323    DO I=1,R22
2324       s11=s1*s11
2325    ENDDO
2326
2327    IF(R2.LT.0) THEN
2328       ! if(old) then
2329       CALL etinv(S11%v%i,S11%v%i)
2330       !       else
2331       !          CALL newetinv(S11%v%j,S11%v%j)
2332       !       endif
2333    ENDIF
2334
2335    powmap=s11
2336
2337
2338    ! powmap=junk
2339    call kill(s11)
2340
2341    master=localmaster
2342
2343  END FUNCTION POWMAP
2344
2345  FUNCTION gPOWMAP( S1, R2 )
2346    implicit none
2347    TYPE (gmap) gPOWMAP
2348    TYPE (gmap), INTENT (IN) :: S1
2349    INTEGER, INTENT (IN) :: R2
2350    TYPE (gmap) S11
2351    INTEGER I,R22
2352    integer localmaster
2353    IF(.NOT.C_%STABLE_DA) RETURN
2354    localmaster=master
2355
2356
2357    gPOWMAP%n=s1%n
2358    call assdamap(gPOWMAP)
2359
2360    call alloc(s11,s1%n)
2361
2362    s11=1
2363
2364
2365    R22=IABS(R2)
2366    DO I=1,R22
2367       s11=s11*s1
2368    ENDDO
2369
2370    IF(R2.LT.0) THEN
2371       ! if(old) then
2372       CALL getinv(S11%v%i,S11%v%i,s11%n)
2373       !       else
2374       !          CALL newgetinv(S11%v%j,S11%v%j,s11%n)
2375       !       endif
2376    ENDIF
2377
2378    gpowmap=s11
2379
2380
2381    ! powmap=junk
2382    call kill(s11)
2383
2384    master=localmaster
2385
2386  END FUNCTION gPOWMAP
2387
2388  FUNCTION gPOWMAPtpsa( S1, R2 )
2389    implicit none
2390    TYPE (gmap) gPOWMAPtpsa
2391    TYPE (gmap), INTENT (IN) :: S1
2392    INTEGER, INTENT (IN) :: R2
2393    TYPE (gmap) S11,s0
2394    INTEGER I,R22
2395    integer localmaster
2396    real(dp), allocatable :: v(:)
2397
2398    IF(.NOT.C_%STABLE_DA) RETURN
2399    localmaster=master
2400
2401
2402    gPOWMAPtpsa%n=s1%n
2403    call assdamap(gPOWMAPtpsa)
2404
2405    call alloc(s11,s1%n)
2406    call alloc(s0,s1%n)
2407
2408    s11=1
2409
2410
2411    R22=IABS(R2)
2412    DO I=1,R22
2413       s11=s11.o.s1
2414    ENDDO
2415    allocate(v(s1%n))
2416
2417    DO I=1,s1%n
2418       v(i)=s11%v(i).sub.'0'
2419    ENDDO
2420
2421
2422    IF(R2.LT.0) THEN
2423       ! if(old) then
2424       CALL getinv(S11%v%i,S11%v%i,s11%n)
2425       !       else
2426       !          CALL newgetinv(S11%v%j,S11%v%j,s11%n)
2427       !       endif
2428    ENDIF
2429
2430    do i=1,s1%n
2431       s0%v(i)=(1.0_dp.mono.i)-v(i)
2432       !s11%v(i)=s11%v(i)-(s11%v(i).sub.'0')
2433    enddo
2434
2435    !     s0=s11.o.s0
2436    s11=s11.o.s0
2437
2438    do i=1,s1%n
2439       gPOWMAPtpsa%v(i)=s11%v(i) !+ (s0%v(i).sub.'0')
2440    enddo
2441
2442
2443
2444
2445
2446
2447    ! powmap=junk
2448    call kill(s11,s0)
2449    deallocate(v)
2450
2451    master=localmaster
2452
2453  END FUNCTION gPOWMAPtpsa
2454
2455  FUNCTION POWMAP_INV( S1, R2 )
2456    implicit none
2457    TYPE (damap) POWMAP_INV
2458    TYPE (damap), INTENT (IN) :: S1
2459    INTEGER, INTENT (IN) :: R2(:)
2460    TYPE (damap) S11
2461    INTEGER I,jn(lnv)
2462    integer localmaster
2463    IF(.NOT.C_%STABLE_DA) RETURN
2464    localmaster=master
2465
2466
2467    do i=1,lnv
2468       jn(i)=0
2469    enddo
2470    do i=1,nd2
2471       jn(i)=R2(I)
2472    enddo
2473    call checkdamap(s1)
2474
2475    call assdamap(POWMAP_INV)
2476
2477    call alloc(s11)
2478
2479    ! if(old) then
2480    if(s1%v(1)%i==0) call crap1("POWMAP_INV 2")  !call etall(s2%m%v%i,nd2)
2481    call etpin(S1%V%i,S11%v%i,jn)
2482    !    else
2483    !       if(.NOT. ASSOCIATED(s1%v(1)%J%r)) call crap1("POWMAP_INV 4")  !call newetall(s2%m%v%J,nd2)
2484    !       call newetpin(S1%v%j,s11%v%j,jn)
2485    !   endif
2486
2487    POWMAP_INV=s11
2488
2489
2490    ! powmap=junk
2491    call kill(s11)
2492
2493    master=localmaster
2494
2495  END FUNCTION POWMAP_INV
2496
2497  subroutine checksymp(s1,norm,orthogonal)
2498    implicit none
2499    TYPE (damap) s1
2500    real(dp)  norm1,mat(8,8),xj(8,8)
2501    real(dp), optional :: norm
2502    integer i,j
2503    logical(lp), optional :: orthogonal
2504    logical(lp) nn
2505    ! checks symplectic conditions on linear map
2506    nn=.not.present(norm)
2507    mat=0.d0
2508    mat=s1
2509    xj=0.d0
2510    if(present(orthogonal)) then
2511       if(orthogonal) then
2512          do i=1,nd
2513             xj(2*i-1,2*i-1)=1.0_dp
2514             xj(2*i,2*i)=1.0_dp
2515          enddo
2516       else
2517          do i=1,nd
2518             xj(2*i-1,2*i)=1.0_dp
2519             xj(2*i,2*i-1)=-1.0_dp
2520          enddo
2521       endif
2522    else
2523       do i=1,nd
2524          xj(2*i-1,2*i)=1.0_dp
2525          xj(2*i,2*i-1)=-1.0_dp
2526       enddo
2527    endif
2528    xj= MATMUL( transpose(mat),MATMUL(xj,mat))
2529
2530    norm1=0.d0
2531    do i=1,nd2
2532       if(lielib_print(9)==1.or.nn)     write(6,'(6(1x,E15.8))') xj(i,1:nd2)
2533       do j=1,nd2
2534
2535          norm1=norm1+abs(xj(i,j))
2536       enddo
2537    enddo
2538    norm1=norm1-nd2
2539    if(lielib_print(9)==1.or.nn) write(6,'(a29,(1x,E15.8))')"deviation from symplecticity ", norm1
2540    if(present(norm)) norm=norm1
2541
2542  end subroutine checksymp
2543
2544
2545  subroutine checkmap(s1)
2546    implicit none
2547    TYPE (damap) s1
2548    integer i
2549    ! if(old) then
2550    do i=1,nd2
2551       if(s1%v(i)%i==0) then
2552          w_p=0
2553          w_p%nc=1
2554          w_p=(/"Should not be here: Assign variables"/)
2555          w_p%fc='(1((1X,A72),/))'
2556          ! call !write_e(200)
2557          !             s1%v(i)%i=dummymap(i)
2558       endif
2559    enddo
2560    !    else
2561    !       do i=1,nd2
2562    !          if(.NOT. ASSOCIATED(s1%v(i)%J%r)) then
2563    !             !             s1%v(i)%j=dummymapl(i)
2564    !             w_p=0
2565    !             w_p%nc=1
2566    !             w_p=(/"Should not be here: Assign variables"/)
2567    !             w_p%fc='(1((1X,A72),/))'
2568    !             ! call !write_e(201)
2569    !         endif
2570    !      enddo
2571    !   endif
2572
2573  end subroutine checkmap
2574
2575  subroutine checkvec(s1)
2576    implicit none
2577    TYPE (vecfield) s1
2578    integer i
2579
2580    ! if(old) then
2581    do i=1,nd2
2582       if(s1%v(i)%i==0) then
2583          !             s1%v(i)%i=dummymap(i)
2584          w_p=0
2585          w_p%nc=1
2586          w_p=(/"Should not be here: Assign variables"/)
2587          w_p%fc='(1((1X,A72),/))'
2588          ! call !write_e(202)
2589       endif
2590    enddo
2591    !    else
2592    !       do i=1,nd2
2593    !          if(.NOT. ASSOCIATED(s1%v(i)%J%r)) then
2594    !             w_p=0
2595    !             w_p%nc=1
2596    !             w_p=(/"Should not be here: Assign variables"/)
2597    !             w_p%fc='(1((1X,A72),/))'
2598    !             ! call !write_e(203)
2599    !             !             s1%v(i)%j=dummymapl(i)
2600    !          endif
2601    !       enddo
2602    !    endif
2603
2604  end subroutine checkvec
2605
2606  subroutine checkpb(s1)
2607    implicit none
2608    TYPE (pbfield) s1
2609
2610    ! if(old) then
2611    if(s1%h%i==0) then
2612       w_p=0
2613       w_p%nc=1
2614       w_p=(/"Should not be here: Assign variables"/)
2615       w_p%fc='(1((1X,A72),/))'
2616       ! call !write_e(204)
2617       !          s1%h%i=dummy
2618    endif
2619    !    else
2620    !       if(.NOT. ASSOCIATED(s1%h%J%r)) then
2621    !          !          s1%h%j=dummyl
2622    !          w_p=0
2623    !          w_p%nc=1
2624    !          w_p=(/"Should not be here: Assign variables"/)
2625    !          w_p%fc='(1((1X,A72),/))'
2626    !          ! call !write_e(205)
2627    !       endif
2628    !   endif
2629  end subroutine checkpb
2630
2631  subroutine checktaylor(s1)
2632    implicit none
2633    TYPE (taylor) s1
2634    ! if(old) then
2635    if(s1%i==0) then
2636       w_p=0
2637       w_p%nc=1
2638       w_p=(/"Should not be here: Assign variables"/)
2639       w_p%fc='(1((1X,A72),/))'
2640       ! call !write_e(206)
2641       !          s1%i=dummy
2642    endif
2643    !    else
2644    !       if(.NOT. ASSOCIATED(s1%J%r)) then
2645    !          w_p=0
2646    !          w_p%nc=1
2647    !          w_p=(/"Should not be here: Assign variables"/)
2648    !          w_p%fc='(1((1X,A72),/))'
2649    !          ! call !write_e(207)
2650    !          !          s1%j=dummyl
2651    !       endif
2652    !   endif
2653  end subroutine checktaylor
2654
2655
2656
2657  subroutine assPB(s1)
2658    implicit none
2659    TYPE (PBFIELD) s1
2660
2661    select case(master)
2662    case(0:ndumt-1)
2663       master=master+1
2664    case(ndumt)
2665       w_p=0
2666       w_p%nc=1
2667       w_p=(/" cannot indent anymore "/)
2668       w_p%fc='(1((1X,A72),/))'
2669       ! call !write_e(100)
2670    end select
2671
2672    call ass0(s1%h)
2673
2674  end subroutine assPB
2675
2676
2677  subroutine asstaylor(s1)
2678    implicit none
2679    TYPE (taylor) s1
2680
2681    select case(master)
2682    case(0:ndumt-1)
2683       master=master+1
2684    case(ndumt)
2685       w_p=0
2686       w_p%nc=1
2687       w_p=(/" cannot indent anymore "/)
2688       w_p%fc='(1((1X,A72),/))'
2689       ! call !write_e(100)
2690    end select
2691
2692    call ass0(s1)
2693
2694  end subroutine asstaylor
2695
2696  subroutine assvec(s1)
2697    implicit none
2698    TYPE (vecfield) s1
2699    integer i
2700
2701    select case(master)
2702    case(0:ndumt-1)
2703       master=master+1
2704    case(ndumt)
2705       w_p=0
2706       w_p%nc=1
2707       w_p=(/" cannot indent anymore "/)
2708       w_p%fc='(1((1X,A72),/))'
2709       ! call !write_e(100)
2710    end select
2711
2712    do i=1,nd2
2713       call ass0(s1%v(i))
2714    enddo
2715
2716  end subroutine assvec
2717
2718  subroutine assmap(s1)
2719    implicit none
2720    TYPE (damap) s1
2721    integer i
2722
2723    select case(master)
2724    case(0:ndumt-1)
2725       master=master+1
2726    case(ndumt)
2727       w_p=0
2728       w_p%nc=1
2729       w_p=(/" cannot indent anymore "/)
2730       w_p%fc='(1((1X,A72),/))'
2731       ! call !write_e(100)
2732    end select
2733
2734    do i=1,nd2
2735       call ass0(s1%v(i))
2736    enddo
2737
2738  end subroutine assmap
2739
2740  subroutine assgmap(s1)
2741    implicit none
2742    TYPE (gmap) s1
2743    integer i
2744
2745    select case(master)
2746    case(0:ndumt-1)
2747       master=master+1
2748    case(ndumt)
2749       w_p=0
2750       w_p%nc=1
2751       w_p=(/" cannot indent anymore "/)
2752       w_p%fc='(1((1X,A72),/))'
2753       ! call !write_e(100)
2754    end select
2755
2756    do i=1,s1%n
2757       call ass0(s1%v(i))
2758    enddo
2759
2760  end subroutine assgmap
2761
2762  !Radiation
2763  SUBROUTINE  allocrad(S1)
2764    implicit none
2765    type (radtaylor),INTENT(INOUT)::S1
2766    ! if(old) then
2767    call etall1(s1%v%i)
2768    call etall(s1%e%i,nd2)
2769    !    else
2770    !       call newetall(s1%v%j,1)
2771    !       call newetall(s1%e%j,nd2)
2772    !   endif
2773  END SUBROUTINE allocrad
2774
2775  SUBROUTINE  KILLrad(S1)
2776    implicit none
2777    type (radtaylor),INTENT(INOUT)::S1
2778    ! if(old) then
2779    call dadal1(s1%v%i)
2780    call dadal(s1%e%i,nd2)
2781    !    else
2782    !       call newdadal(s1%v%j,1)
2783    !       call newdadal(s1%e%j,nd2)
2784    !   endif
2785
2786  END SUBROUTINE KILLrad
2787
2788  SUBROUTINE  allocrads(S1,n)
2789    implicit none
2790    type (radtaylor),INTENT(INOUT),dimension(:)::S1
2791    integer i ,n
2792    do i=1,n
2793       call allocrad(s1(i))
2794    enddo
2795
2796  END SUBROUTINE allocrads
2797
2798  SUBROUTINE  KILLrads(S1,n)
2799    implicit none
2800    type (radtaylor),INTENT(INOUT),dimension(:)::s1
2801    integer i,n
2802    do i=1,n
2803       call killrad(s1(i))
2804    enddo
2805
2806
2807  END SUBROUTINE KILLrads
2808
2809
2810
2811  SUBROUTINE  EQUALrad(S2,S1)
2812    implicit none
2813    type (TAYLOR),INTENT(inOUT)::S2
2814    type (radTAYLOR),INTENT(IN)::S1
2815    !     if(iass0>ndum) then
2816    !     call ndum_warning
2817    !     endif
2818    !      iass0=0
2819    ! if(old) then
2820    if(s2%i==0) call crap1("EQUALrad 1")  ! call allocw_old(s2)
2821    CALL DACOP(S1%v%I,S2%I)
2822    !    else
2823    !       IF (.NOT. ASSOCIATED(s2%j%r)) call crap1("EQUALrad 2")  !call allocw_old(s2)
2824    !       call newdacop(S1%v%j,S2%j)
2825    !   endif
2826  END SUBROUTINE EQUALrad
2827
2828
2829  SUBROUTINE  radEQUAL(S2,S1)
2830    implicit none
2831    type (radTAYLOR),INTENT(inOUT)::S2
2832    type (TAYLOR),INTENT(IN)::S1
2833    call check_snake
2834
2835    !     if(iass0>ndum) then
2836    !     call ndum_warning
2837    !     endif
2838    !      iass0=0
2839    ! if(old) then
2840    if(s2%v%i==0) call crap1("radEQUAL 1")  ! call allocw_old(s2%v)
2841    CALL DACOP(S1%I,S2%v%I)
2842    !    else
2843    !       IF (.NOT. ASSOCIATED(s2%v%j%r)) call crap1("radEQUAL 2")  !call allocw_old(s2%v)
2844    !       call newdacop(S1%j,S2%v%j)
2845    !   endif
2846  END SUBROUTINE radEQUAL
2847
2848  ! remove small numbers
2849
2850
2851  SUBROUTINE  flip_damap(S1,S2)
2852    implicit none
2853    type (damap),INTENT(INOUT)::S2
2854    type (damap), intent(INOUT):: s1
2855
2856    call flip(S1%v%i,S2%v%i)
2857
2858  end SUBROUTINE  flip_damap
2859
2860  SUBROUTINE  flip_taylor(S1,S2,i)
2861    implicit none
2862    type (taylor),INTENT(INOUT)::S2
2863    type (taylor), intent(INOUT):: s1
2864    integer i
2865    call flip_i(S1%i,S2%i,i)
2866
2867  end SUBROUTINE  flip_taylor
2868
2869  SUBROUTINE  flip_vecfield(S1,S2,i)
2870    implicit none
2871    type (vecfield),INTENT(INOUT)::S2
2872    type (vecfield), intent(INOUT):: s1
2873    integer i
2874    call flipflo(S1%v%i,S2%v%i,i)
2875
2876  end SUBROUTINE  flip_vecfield
2877
2878end module tpsalie
Note: See TracBrowser for help on using the repository browser.