1      subroutine tce_print_x1(d_r1,k_r1_offset,r1_thresh,irrep)
2c
3c $Id$
4c
5c Prints the X1 tensor
6c
7      implicit none
8#include "global.fh"
9#include "mafdecls.fh"
10#include "sym.fh"
11#include "util.fh"
12#include "stdio.fh"
13#include "errquit.fh"
14#include "tce.fh"
15#include "tce_main.fh"
16      integer d_r1
17      integer p1b
18      integer h2b
19      integer p1
20      integer h2
21      integer k_r1_offset
22      integer size
23      integer l_r1,k_r1
24      integer i
25      integer irrep
26      logical nodezero
27      double precision r1_thresh
28      character*4 irrep_p,irrep_h
29      character*5 spinnames(2),spin_p,spin_h
30      data spinnames/'alpha','beta '/
31c
32      nodezero = (ga_nodeid().eq.0)
33      if (nodezero) then
34         write(LuOut,9000) 'Singles contributions'
35      do p1b = noab+1,noab+nvab
36        do h2b = 1,noab
37          if (int_mb(k_spin+p1b-1) .eq. int_mb(k_spin+h2b-1)) then
38            if ((.not.restricted).or.(int_mb(k_spin+p1b-1)
39     1        +int_mb(k_spin+h2b-1).ne.4)) then
40              if (ieor(int_mb(k_sym+p1b-1),int_mb(k_sym+h2b-1))
41     1          .eq.irrep) then
42                size = int_mb(k_range+p1b-1) * int_mb(k_range+h2b-1)
43                if (.not.ma_push_get(mt_dbl,size,'r1',l_r1,k_r1))
44     1            call errquit('tce_print_x1: MA problem',0,MA_ERR)
45                call get_hash_block(d_r1,dbl_mb(k_r1),size,
46     1            int_mb(k_r1_offset),((p1b-noab-1)*noab+h2b-1))
47                i = 0
48                do p1 = 1,int_mb(k_range+p1b-1)
49                  do h2 = 1,int_mb(k_range+h2b-1)
50                    i = i + 1
51                    if (dabs(dbl_mb(k_r1+i-1)).gt.r1_thresh) then
52                      call sym_irrepname(geom,int_mb(k_sym+p1b-1)+1,
53     1                                   irrep_p)
54                      call sym_irrepname(geom,int_mb(k_sym+h2b-1)+1,
55     1                                   irrep_h)
56                      spin_p = spinnames(int_mb(k_spin+p1b-1))
57                      spin_h = spinnames(int_mb(k_spin+h2b-1))
58                      if (nodezero) write(LuOut,9010)
59c     1                int_mb(k_offset+p1b-1)+p1,irrep_p,spin_p,
60c     2                int_mb(k_offset+h2b-1)+h2,irrep_h,spin_h,
61c     3                dbl_mb(k_r1+i-1)
62     1                (int_mb(k_mo_index+int_mb(k_offset+p1b-1)+p1-1)+
63     1                 2-int_mb(k_spin+p1b-1))/2,
64     1                irrep_p,spin_p,
65     2                (int_mb(k_mo_index+int_mb(k_offset+h2b-1)+h2-1)+
66     1                 2-int_mb(k_spin+h2b-1))/2,
67     2                irrep_h,spin_h,
68     2                dbl_mb(k_r1+i-1)
69                    endif
70                  enddo
71                enddo
72                if (.not.ma_pop_stack(l_r1))
73     1            call errquit('tce_print_x1: MA problem',1,MA_ERR)
74              endif
75            endif
76          endif
77        enddo
78      enddo
79      endif
80      call ga_sync()
81 9000 format(/,1x,a)
82 9010 format(1x,i5,a4,'(',a5,') --- ',i5,a4,'(',a5,')',f20.10)
83      return
84      end
85
86      subroutine tce_print_y1(d_r1,k_r1_offset,r1_thresh,irrep)
87c
88c Prints the Y1 tensor
89c
90      implicit none
91#include "global.fh"
92#include "mafdecls.fh"
93#include "sym.fh"
94#include "util.fh"
95#include "stdio.fh"
96#include "errquit.fh"
97#include "tce.fh"
98#include "tce_main.fh"
99      integer d_r1
100      integer p1b
101      integer h2b
102      integer p1
103      integer h2
104      integer k_r1_offset
105      integer size
106      integer l_r1,k_r1
107      integer i
108      integer irrep
109      logical nodezero
110      double precision r1_thresh
111      character*4 irrep_p,irrep_h
112      character*5 spinnames(2),spin_p,spin_h
113      data spinnames/'alpha','beta '/
114c
115      nodezero = (ga_nodeid().eq.0)
116      if (nodezero) then
117         write(LuOut,9000) 'Singles contributions'
118      do h2b = 1,noab
119        do p1b = noab+1,noab+nvab
120          if (int_mb(k_spin+p1b-1) .eq. int_mb(k_spin+h2b-1)) then
121          if ((.not.restricted).or.(int_mb(k_spin+p1b-1)
122     1      +int_mb(k_spin+h2b-1).ne.4)) then
123          if (ieor(int_mb(k_sym+p1b-1),int_mb(k_sym+h2b-1))
124     1      .eq.irrep) then
125            size = int_mb(k_range+p1b-1) * int_mb(k_range+h2b-1)
126            if (.not.ma_push_get(mt_dbl,size,'r1',l_r1,k_r1))
127     1        call errquit('tce_print_y1: MA problem',0,MA_ERR)
128            call get_hash_block(d_r1,dbl_mb(k_r1),size,
129     1        int_mb(k_r1_offset),((h2b-1)*nvab+p1b-noab-1))
130            i = 0
131            do p1 = 1,int_mb(k_range+p1b-1)
132              do h2 = 1,int_mb(k_range+h2b-1)
133                i = i + 1
134                if (dabs(dbl_mb(k_r1+i-1)).gt.r1_thresh) then
135                  call sym_irrepname(geom,int_mb(k_sym+p1b-1)+1,irrep_p)
136                  call sym_irrepname(geom,int_mb(k_sym+h2b-1)+1,irrep_h)
137                  spin_p = spinnames(int_mb(k_spin+p1b-1))
138                  spin_h = spinnames(int_mb(k_spin+h2b-1))
139                  if (nodezero) write(LuOut,9010)
140     1              int_mb(k_offset+p1b-1)+p1,irrep_p,spin_p,
141     2              int_mb(k_offset+h2b-1)+h2,irrep_h,spin_h,
142     3              dbl_mb(k_r1+i-1)
143                endif
144              enddo
145            enddo
146            if (.not.ma_pop_stack(l_r1))
147     1        call errquit('tce_print_y1: MA problem',1,MA_ERR)
148          endif
149          endif
150          endif
151        enddo
152      enddo
153      endif
154      call ga_sync()
155 9000 format(/,1x,a)
156 9010 format(1x,i5,a4,'(',a5,') --- ',i5,a4,'(',a5,')',f20.10)
157      return
158      end
159c
160c
161c
162c
163c
164      subroutine tce_print_g1(rtdb,d_r1,k_r1_offset,r1_thresh,irrep_r1)
165c
166c Prints the G1 tensor
167c
168      implicit none
169#include "global.fh"
170#include "mafdecls.fh"
171#include "sym.fh"
172#include "util.fh"
173#include "stdio.fh"
174#include "errquit.fh"
175#include "rtdb.fh"
176#include "tce.fh"
177#include "tce_main.fh"
178      integer rtdb
179      integer d_r1
180      integer g1b
181      integer g2b
182      integer g1
183      integer g2
184      integer k_r1_offset
185      integer size
186      integer l_r1,k_r1
187      integer i
188      integer irrep_r1
189      integer eiformat
190      logical nodezero
191      double precision r1_thresh
192      character*4 irrep_g1,irrep_g2
193      character*5 spinnames(2),spin_g1,spin_g2
194      data spinnames/'alpha','beta '/
195c
196      if (.not.rtdb_get(rtdb,'tce:eiformat',mt_int,1,eiformat)) then
197        eiformat = 1
198        if (.not.rtdb_put(rtdb,'tce:eiformat',mt_int,1,eiformat))
199     1  call errquit('tce_init: failed to write eiformat',0,RTDB_ERR)
200      endif
201c
202      nodezero = (ga_nodeid().eq.0)
203c      if (nodezero) write(LuOut,9000) 'Singles contributions'
204      if(nodezero) then
205      do g2b = 1,noab+nvab
206        do g1b = 1,noab+nvab
207          if (int_mb(k_spin+g1b-1) .eq. int_mb(k_spin+g2b-1)) then
208          if ((.not.restricted).or.(int_mb(k_spin+g1b-1)
209     1      +int_mb(k_spin+g2b-1).ne.4)) then
210          if (ieor(int_mb(k_sym+g1b-1),int_mb(k_sym+g2b-1))
211     1      .eq.irrep_r1) then
212            size = int_mb(k_range+g1b-1) * int_mb(k_range+g2b-1)
213            if (.not.ma_push_get(mt_dbl,size,'r1',l_r1,k_r1))
214     1        call errquit('tce_print_y1: MA problem',0,MA_ERR)
215            call get_hash_block(d_r1,dbl_mb(k_r1),size,
216     1        int_mb(k_r1_offset),(g2b-1)*(noab+nvab)+g1b-1)
217            i = 0
218            do g1 = 1,int_mb(k_range+g1b-1)
219             do g2 = 1,int_mb(k_range+g2b-1)
220               i = i + 1
221               if (dabs(dbl_mb(k_r1+i-1)).gt.r1_thresh) then
222                 call sym_irrepname(geom,int_mb(k_sym+g1b-1)+1,irrep_g1)
223                 call sym_irrepname(geom,int_mb(k_sym+g2b-1)+1,irrep_g2)
224                 spin_g1 = spinnames(int_mb(k_spin+g1b-1))
225                 spin_g2 = spinnames(int_mb(k_spin+g2b-1))
226                 if ((nodezero).and.(eiformat.eq.1)) then
227                   write(LuOut,9011)
228     1             int_mb(k_offset+g1b-1)+g1,irrep_g1,spin_g1,
229     2             int_mb(k_offset+g2b-1)+g2,irrep_g2,spin_g2,
230     3             dbl_mb(k_r1+i-1)
231                 elseif ((nodezero).and.(eiformat.eq.2)) then
232                   write(LuOut,9012)
233     1             int_mb(k_offset+g1b-1)+g1,irrep_g1,spin_g1,
234     2             int_mb(k_offset+g2b-1)+g2,irrep_g2,spin_g2,
235     3             dbl_mb(k_r1+i-1)
236                 elseif ((nodezero).and.(eiformat.eq.3)) then
237                   write(LuOut,9013)
238     1             int_mb(k_offset+g1b-1)+g1,int_mb(k_sym+g1b-1)+1,
239     2             int_mb(k_spin+g1b-1),
240     3             int_mb(k_offset+g2b-1)+g2,int_mb(k_sym+g2b-1)+1,
241     4             int_mb(k_spin+g2b-1),
242     5             dbl_mb(k_r1+i-1)
243                 elseif ((nodezero).and.(eiformat.eq.4)) then
244                   write(LuOut,9014)
245     1             int_mb(k_offset+g1b-1)+g1,
246     2             int_mb(k_offset+g2b-1)+g2,
247     3             dbl_mb(k_r1+i-1)
248                 endif
249               endif
250             enddo
251            enddo
252            if (.not.ma_pop_stack(l_r1))
253     1        call errquit('tce_print_y1: MA problem',1,MA_ERR)
254          endif
255          endif
256          endif
257        enddo
258      enddo
259      endif
260      call ga_sync()
261 9000 format(/,1x,a)
262 9011 format(1x,i3,a4,'(',a5,') --- ',i3,a4,'(',a5,')',f20.14)
263 9012 format(1x,i3,',  ',a4,',',a5,',',i3,',  ',a4,',',a5,',',f20.14)
264 9013 format(1x,i3,',',i3,',',i3,',',i3,',',i3,',',i3,',',f20.14)
265 9014 format(1x,i3,', ',i3,', ',f20.14)
266      return
267      end
268c
269      subroutine tce_string_x1(d_r1,k_r1_offset,r1_thresh,irrep,xnorm)
270c
271c $Id: tce_print_x1.F 26707 2015-01-17 01:06:43Z edo $
272c
273c Prints the X1 tensor
274c
275      implicit none
276#include "global.fh"
277#include "mafdecls.fh"
278#include "sym.fh"
279#include "util.fh"
280#include "stdio.fh"
281#include "errquit.fh"
282#include "tce.fh"
283#include "tce_main.fh"
284      integer d_r1
285      integer p1b
286      integer h2b
287      integer p1
288      integer h2
289      integer k_r1_offset
290      integer size
291      integer l_r1,k_r1
292      integer i
293      integer irrep
294c string
295       integer spinsum,j,iposh,iposp
296       integer nosa,nosb,nvsa,nvsb
297c ------
298      logical nodezero
299      double precision r1_thresh,xnorm
300      character*4 irrep_p,irrep_h
301      character*5 spinnames(2),spin_p,spin_h
302c string
303      character*6 refstring(202)
304      character*1 bra,num1
305      character*2 num2
306      character*3 keta,ketb,kea,keb
307      character*4 cre4
308      character*5 cre5
309      character*6 cre6
310c ------
311      data spinnames/'alpha','beta '/
312c
313c
314      nodezero = (ga_nodeid().eq.0)
315      if (nodezero) then
316cc        write(6,*)'check 2: xnorm =',xnorm
317cc        call util_flush(6)
318c
319c string
320c
321        bra='('
322        keta='a)+'
323        ketb='b)+'
324        kea='a)'
325        keb='b)'
326c
327        nosa=0
328        nosb=0
329        nvsa=0
330        nvsb=0
331        do i=1,noa
332         nosa=nosa+int_mb(k_range+i-1)
333        enddo
334        do i=1,nob
335         nosb=nosb+int_mb(k_range+noa+i-1)
336        enddo
337        do i=1,nva
338         nvsa=nvsa+int_mb(k_range+noab+i-1)
339        enddo
340        do i=1,nvb
341         nvsb=nvsb+int_mb(k_range+noab+nva+i-1)
342        enddo
343c reference string
344        do i=1,nosa
345          if(i.le.9) then
346           j=i
347           write(num1,301) j
348           cre5=bra//num1//keta
349           refstring(i+2)=cre5
350          else
351           j=i
352           write(num2,302) j
353           cre6=bra//num2//keta
354           refstring(i+2)=cre6
355          endif
356        enddo
357c
358        do i=1,nosb
359          if(i.le.9) then
360           j=i
361           write(num1,301) j
362           cre5=bra//num1//ketb
363           refstring(i+nosa+2)=cre5
364          else
365           j=i
366           write(num2,302) j
367           cre6=bra//num2//ketb
368           refstring(i+nosa+2)=cre6
369          endif
370        enddo
371c
372c ------
373         write(LuOut,9000) 'Singles strings'
374      do p1b = noab+1,noab+nvab
375        do h2b = 1,noab
376          if (int_mb(k_spin+p1b-1) .eq. int_mb(k_spin+h2b-1)) then
377            spinsum=int_mb(k_spin+p1b-1)+int_mb(k_spin+h2b-1)
378            if ((.not.restricted).or.(int_mb(k_spin+p1b-1)
379     1        +int_mb(k_spin+h2b-1).ne.4)) then
380              if (ieor(int_mb(k_sym+p1b-1),int_mb(k_sym+h2b-1))
381     1          .eq.irrep) then
382                size = int_mb(k_range+p1b-1) * int_mb(k_range+h2b-1)
383                if (.not.ma_push_get(mt_dbl,size,'r1',l_r1,k_r1))
384     1            call errquit('tce_print_x1: MA problem',0,MA_ERR)
385                call get_hash_block(d_r1,dbl_mb(k_r1),size,
386     1            int_mb(k_r1_offset),((p1b-noab-1)*noab+h2b-1))
387                i = 0
388                do p1 = 1,int_mb(k_range+p1b-1)
389                  do h2 = 1,int_mb(k_range+h2b-1)
390                    i = i + 1
391                    if (dabs(dbl_mb(k_r1+i-1)).gt.r1_thresh) then
392                     iposh=int_mb(k_offset+h2b-1)+h2
393                     iposp=int_mb(k_offset+p1b-1)+p1-nosb
394c
395                 if(iposh.le.9) then
396                  write(num1,301) iposh
397                  cre4=bra//num1//kea
398                  refstring(2)=cre4
399                 else
400                  write(num2,302) iposh
401                  cre5=bra//num2//kea
402                  refstring(2)=cre5
403                 endif
404c
405                 if(iposp.le.9) then
406                  write(num1,301) iposp
407                  cre5=bra//num1//keta
408                  refstring(1)=cre5
409                 else
410                  write(num2,302) iposp
411                  cre6=bra//num2//keta
412                  refstring(1)=cre6
413                 endif
414c yaml initial_state_suggestion state
415c yaml "Singles contribution"
416        write(6,*) xnorm*dbl_mb(k_r1+i-1),':',
417     &             (refstring(j),j=1,nosa+nosb+2),
418     &             '|0>'
419c spinsum.eq.4 - mirroring beta->beta excitation
420             if(spinsum.eq.2) then
421                 if(iposh.le.9) then
422                  write(num1,301) iposh
423                  cre4=bra//num1//keb
424                  refstring(2)=cre4
425                 else
426                  write(num2,302) iposh
427                  cre5=bra//num2//keb
428                  refstring(2)=cre5
429                 endif
430c
431                 if(iposp.le.9) then
432                  write(num1,301) iposp
433                  cre5=bra//num1//ketb
434                  refstring(1)=cre5
435                 else
436                  write(num2,302) iposp
437                  cre6=bra//num2//ketb
438                  refstring(1)=cre6
439                 endif
440c yaml initial_state_suggestion state
441c yaml "Singles contribution"
442        write(6,*) xnorm*dbl_mb(k_r1+i-1),':',
443     &             (refstring(j),j=1,nosa+nosb+2),
444     &             '|0>'
445             endif  !spinsum
446c ------
447c
448c
449                    endif
450                  enddo
451                enddo
452                if (.not.ma_pop_stack(l_r1))
453     1            call errquit('tce_print_x1: MA problem',1,MA_ERR)
454              endif
455            endif
456          endif
457        enddo
458      enddo
459      endif
460      call ga_sync()
461  301 format(i1)
462  302 format(i2)
463 9000 format(/,1x,a)
464 9010 format(1x,i5,a4,'(',a5,') --- ',i5,a4,'(',a5,')',f20.10)
465      return
466      end
467c
468c
469c
470c
471c
472      subroutine tce_string_x1_norm(d_r1,k_r1_offset,r1_thresh,irrep,
473     &           xnorm1)
474c
475c $Id: tce_print_x1.F 26707 2015-01-17 01:06:43Z edo $
476c
477c Prints the X1 tensor
478c
479      implicit none
480#include "global.fh"
481#include "mafdecls.fh"
482#include "sym.fh"
483#include "util.fh"
484#include "stdio.fh"
485#include "errquit.fh"
486#include "tce.fh"
487#include "tce_main.fh"
488      integer d_r1
489      integer p1b
490      integer h2b
491      integer p1
492      integer h2
493      integer k_r1_offset
494      integer size
495      integer l_r1,k_r1
496      integer i
497      integer irrep
498c string
499       integer spinsum,j
500c ------
501      logical nodezero
502      double precision r1_thresh,xnorm1
503      character*4 irrep_p,irrep_h
504      character*5 spinnames(2),spin_p,spin_h
505c ------
506      data spinnames/'alpha','beta '/
507c xnorm defined only on nodezero
508      xnorm1=0.0d0
509c
510      nodezero = (ga_nodeid().eq.0)
511      if (nodezero) then
512c
513         write(LuOut,9000) 'singles norm'
514      do p1b = noab+1,noab+nvab
515        do h2b = 1,noab
516          if (int_mb(k_spin+p1b-1) .eq. int_mb(k_spin+h2b-1)) then
517            spinsum=int_mb(k_spin+p1b-1)+int_mb(k_spin+h2b-1)
518            if ((.not.restricted).or.(int_mb(k_spin+p1b-1)
519     1        +int_mb(k_spin+h2b-1).ne.4)) then
520              if (ieor(int_mb(k_sym+p1b-1),int_mb(k_sym+h2b-1))
521     1          .eq.irrep) then
522                size = int_mb(k_range+p1b-1) * int_mb(k_range+h2b-1)
523                if (.not.ma_push_get(mt_dbl,size,'r1',l_r1,k_r1))
524     1            call errquit('tce_print_x1: MA problem',0,MA_ERR)
525                call get_hash_block(d_r1,dbl_mb(k_r1),size,
526     1            int_mb(k_r1_offset),((p1b-noab-1)*noab+h2b-1))
527                i = 0
528                do p1 = 1,int_mb(k_range+p1b-1)
529                  do h2 = 1,int_mb(k_range+h2b-1)
530                    i = i + 1
531                    if (dabs(dbl_mb(k_r1+i-1)).gt.r1_thresh) then
532c
533c string for RHF reference only  a->a
534                 xnorm1=xnorm1+dbl_mb(k_r1+i-1)*dbl_mb(k_r1+i-1)
535c spinsum.eq.4 - mirroring beta->beta excitation
536             if(spinsum.eq.2) then
537                 xnorm1=xnorm1+dbl_mb(k_r1+i-1)*dbl_mb(k_r1+i-1)
538             endif  !spinsum
539c ------
540c
541c
542                    endif
543                  enddo
544                enddo
545                if (.not.ma_pop_stack(l_r1))
546     1            call errquit('tce_print_x1: MA problem',1,MA_ERR)
547              endif
548            endif
549          endif
550        enddo
551      enddo
552      endif
553      call ga_sync()
554 9000 format(/,1x,a)
555 9010 format(1x,i5,a4,'(',a5,') --- ',i5,a4,'(',a5,')',f20.10)
556      return
557      end
558c
559c
560c
561c
562c
563c
564      subroutine tce_string_ref(xnorm)
565c
566c $Id: tce_print_x1.F 26707 2015-01-17 01:06:43Z edo $
567c
568c Prints the X1 tensor
569c
570      implicit none
571#include "global.fh"
572#include "mafdecls.fh"
573#include "sym.fh"
574#include "util.fh"
575#include "stdio.fh"
576#include "errquit.fh"
577#include "tce.fh"
578#include "tce_main.fh"
579      integer p1b
580      integer h2b
581      integer p1
582      integer h2
583      integer i
584c string
585       integer spinsum,j,iposh,iposp
586       integer nosa,nosb,nvsa,nvsb
587c ------
588      logical nodezero
589      double precision xnorm
590      character*4 irrep_p,irrep_h
591      character*5 spinnames(2),spin_p,spin_h
592c string
593      character*6 refstring(202)
594      character*1 bra,num1
595      character*2 num2
596      character*3 keta,ketb,kea,keb
597      character*4 cre4
598      character*5 cre5
599      character*6 cre6
600c ------
601      data spinnames/'alpha','beta '/
602c
603c
604      nodezero = (ga_nodeid().eq.0)
605      if (nodezero) then
606c        write(6,*)'check 2: xnorm =',xnorm
607c        call util_flush(6)
608c
609c string
610c
611        bra='('
612        keta='a)+'
613        ketb='b)+'
614        kea='a)'
615        keb='b)'
616c
617        nosa=0
618        nosb=0
619        nvsa=0
620        nvsb=0
621        do i=1,noa
622         nosa=nosa+int_mb(k_range+i-1)
623        enddo
624        do i=1,nob
625         nosb=nosb+int_mb(k_range+noa+i-1)
626        enddo
627        do i=1,nva
628         nvsa=nvsa+int_mb(k_range+noab+i-1)
629        enddo
630        do i=1,nvb
631         nvsb=nvsb+int_mb(k_range+noab+nva+i-1)
632        enddo
633c reference string
634        do i=1,nosa
635          if(i.le.9) then
636           j=i
637           write(num1,301) j
638           cre5=bra//num1//keta
639           refstring(i+2)=cre5
640          else
641           j=i
642           write(num2,302) j
643           cre6=bra//num2//keta
644           refstring(i+2)=cre6
645          endif
646        enddo
647c
648        do i=1,nosb
649          if(i.le.9) then
650           j=i
651           write(num1,301) j
652           cre5=bra//num1//ketb
653           refstring(i+nosa+2)=cre5
654          else
655           j=i
656           write(num2,302) j
657           cre6=bra//num2//ketb
658           refstring(i+nosa+2)=cre6
659          endif
660        enddo
661c
662         write(LuOut,9000) 'Reference string'
663c yaml initial_state_suggestion state
664         write(6,*) xnorm,':',
665     &             (refstring(j),j=3,nosa+nosb+2),
666     &             '|0>'
667      endif !nodezero
668      call ga_sync()
669  301 format(i1)
670  302 format(i2)
671 9000 format(/,1x,a)
672 9010 format(1x,i5,a4,'(',a5,') --- ',i5,a4,'(',a5,')',f20.10)
673      return
674      end
675c
676c
677c
678c
679c
680