1* $Id$
2*
3
4
5*     ***********************************
6*     *			                *
7*     *		c_geodesic_init	        *
8*     *					*
9*     ***********************************
10*
11*     Uses - c_geodesic common block
12*
13
14      subroutine c_geodesic_init()
15      implicit none
16
17#include "bafdecls.fh"
18#include "errquit.fh"
19#include "c_geodesic_common.fh"
20
21*     **** local variables ****
22      integer npack1,neall,nemax,nbrillq
23
24*     **** external functions ****
25      integer  cpsi_ne,cpsi_neq,cpsi_nbrillq,Pneb_w_size
26      integer  cpsi_data_alloc
27      external cpsi_ne,cpsi_neq,cpsi_nbrillq,Pneb_w_size
28      external cpsi_data_alloc
29
30      call Cram_max_npack(npack1)
31      neall  = cpsi_neq(1)+cpsi_neq(2)
32      nemax  = cpsi_ne(1)+cpsi_ne(2)
33      nbrillq = cpsi_nbrillq()
34
35      U_tag  = cpsi_data_alloc(nbrillq,neall,2*npack1)
36      Vt_tag = cpsi_data_alloc(nbrillq,1,2*Pneb_w_size(0,1))
37      S_tag  = cpsi_data_alloc(nbrillq,nemax,1)
38      return
39      end
40
41*     ***********************************
42*     *					*
43*     *		c_geodesic_finalize	*
44*     *					*
45*     ***********************************
46*
47*     Uses - c_geodesic common block
48*
49      subroutine c_geodesic_finalize()
50      implicit none
51
52#include "bafdecls.fh"
53#include "errquit.fh"
54#include "c_geodesic_common.fh"
55
56
57      call cpsi_data_dealloc(U_tag)
58      call cpsi_data_dealloc(Vt_tag)
59      call cpsi_data_dealloc(S_tag)
60      return
61      end
62
63
64
65*     ***********************************
66*     *					*
67*     *		c_geodesic_start        *
68*     *					*
69*     ***********************************
70*
71*     Uses - c_geodesic common block
72*
73
74      subroutine c_geodesic_start(A_tag,max_sigma,dE)
75      implicit none
76      integer A_tag
77      real*8  max_sigma,dE
78
79#include "bafdecls.fh"
80#include "errquit.fh"
81#include "c_geodesic_common.fh"
82
83*     **** local variables ****
84      integer nb,i,nbrillq,neall,npack1
85      integer ashift,ushift,sshift,vshift,vtshift,V_tag
86      real*8 tmp
87
88*     **** external functions ****
89      integer  cpsi_nbrillq,cpsi_neq,Pneb_w_size
90      integer  cpsi_data_push_stack,cpsi_data_get_chnk
91      real*8   c_electron_eorbit
92      external cpsi_nbrillq,cpsi_neq,Pneb_w_size
93      external cpsi_data_push_stack,cpsi_data_get_chnk
94      external c_electron_eorbit
95
96
97      call nwpw_timing_start(10)
98
99      call Cram_max_npack(npack1)
100      nbrillq = cpsi_nbrillq()
101      neall   = cpsi_neq(1)+cpsi_neq(2)
102
103*     **** allocate tmp space ****
104      V_tag = cpsi_data_push_stack(nbrillq,1,2*Pneb_w_size(0,1))
105
106      max_sigma = -1.0d200
107      do nb=1,nbrillq
108         ashift  = cpsi_data_get_chnk(A_tag,nb)
109         ushift  = cpsi_data_get_chnk(U_tag,nb)
110         vshift  = cpsi_data_get_chnk(V_tag,nb)
111         sshift  = cpsi_data_get_chnk(S_tag,nb)
112         vtshift = cpsi_data_get_chnk(Vt_tag,nb)
113
114*        **** HomeGrown SVD ****
115         call Pneb_SVD(0,nb,npack1,
116     >              dbl_mb(ashift),
117     >              dbl_mb(ushift),
118     >              dbl_mb(sshift),
119     >              dbl_mb(vshift) )
120
121*        **** calculate Vt ****
122         call Pneb_w_copy_dagger(0,nb,dbl_mb(vshift),dbl_mb(vtshift))
123
124*        **** find max_sigma ****
125         do i=1,neall
126            tmp=dabs(dbl_mb(sshift))
127            if (tmp.gt.max_sigma) max_sigma=tmp
128            sshift=sshift+1
129         end do
130
131      end do
132      call K1dB_MaxAll(max_sigma)
133
134*     **** calculate 2*<A|H|psi> ****
135      dE = 2.0d0*c_electron_eorbit(A_tag)
136
137*     **** deallocate tmp space ****
138      call cpsi_data_pop_stack(V_tag)
139
140      call nwpw_timing_end(10)
141      return
142      end
143
144
145
146*     *******************************
147*     *		            	    *
148*     *		c_geodesic_get	    *
149*     *			            *
150*     *******************************
151*
152*     Uses - c_geodesic common block
153*
154
155      subroutine c_geodesic_get(t,Yold_tag,Ynew_tag)
156      implicit none
157      real*8     t
158      integer    Yold_tag
159      integer    Ynew_tag
160
161#include "bafdecls.fh"
162#include "errquit.fh"
163#include "c_geodesic_common.fh"
164
165*     **** local variables ****
166      complex*16 zero,one
167      parameter (zero=(0.0d0,0.0d0), one=(1.0d0,0.0d0))
168
169      logical value
170      integer nb,npack1,nemax,nbrillq
171      integer tmp1(2),tmp2(2),tmp3(2),tmpC(2),tmpS(2)
172      integer yoldshift,ynewshift,ushift,sshift,vtshift
173
174
175*     **** external functions ****
176      logical  Pneb_w_push_get,Pneb_w_pop_stack
177      integer  cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
178      external Pneb_w_push_get,Pneb_w_pop_stack
179      external cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
180
181      call nwpw_timing_start(10)
182      call Cram_max_npack(npack1)
183      nemax   = cpsi_ne(1)+cpsi_ne(2)
184      nbrillq = cpsi_nbrillq()
185
186*     **** allocate tmp space ****
187      value =           Pneb_w_push_get(0,1,tmp1)
188      value = value.and.Pneb_w_push_get(0,1,tmp2)
189      value = value.and.Pneb_w_push_get(0,1,tmp3)
190      value = value.and.BA_push_get(mt_dbl,nemax,'tmpC',tmpC(2),tmpC(1))
191      value = value.and.BA_push_get(mt_dbl,nemax,'tmpS',tmpS(2),tmpS(1))
192      if (.not. value)
193     > call errquit('c_geodesic_get: out of stack memory',0,MA_ERR)
194
195      do nb=1,nbrillq
196         yoldshift = cpsi_data_get_chnk(Yold_tag,nb)
197         ynewshift = cpsi_data_get_chnk(Ynew_tag,nb)
198         ushift    = cpsi_data_get_chnk(U_tag,nb)
199         sshift    = cpsi_data_get_chnk(S_tag,nb)
200         vtshift   = cpsi_data_get_chnk(Vt_tag,nb)
201         call Pneb_SCVtrans1(0,nb,t,
202     >                    dbl_mb(sshift),
203     >                    dbl_mb(vtshift),
204     >                    dcpl_mb(tmp1(1)),
205     >                    dcpl_mb(tmp3(1)),
206     >                    dbl_mb(tmpC(1)),
207     >                    dbl_mb(tmpS(1)))
208
209         call Pneb_www_Multiply2(0,nb,
210     >                    one,
211     >                    dbl_mb(vtshift),
212     >                    dcpl_mb(tmp1(1)),
213     >                    zero,
214     >                    dcpl_mb(tmp2(1)))
215
216         call Pneb_fwf_Multiply(0,nb,
217     >                    one,
218     >                    dbl_mb(yoldshift),npack1,
219     >                    dcpl_mb(tmp2(1)),
220     >                    zero,
221     >                    dbl_mb(ynewshift))
222
223         call Pneb_fwf_Multiply(0,nb,
224     >                    one,
225     >                    dbl_mb(ushift),npack1,
226     >                    dcpl_mb(tmp3(1)),
227     >                    one,
228     >                    dbl_mb(ynewshift))
229
230*        **** Orthonormality Check ****
231         call Pneb_orthoCheckMake(.true.,0,nb,npack1,dbl_mb(ynewshift))
232      end do
233
234*     **** deallocate tmp space ****
235      value =           BA_pop_stack(tmpS(2))
236      value = value.and.BA_pop_stack(tmpC(2))
237      value = value.and.Pneb_w_pop_stack(tmp3)
238      value = value.and.Pneb_w_pop_stack(tmp2)
239      value = value.and.Pneb_w_pop_stack(tmp1)
240      if (.not. value)
241     > call errquit('error popping stack memory',0, MA_ERR)
242
243      call nwpw_timing_end(10)
244
245      return
246      end
247
248*     ***********************************
249*     *					*
250*     *		c_geodesic_transport	*
251*     *					*
252*     ***********************************
253*
254*     Uses - geodesic common block
255*
256
257      subroutine c_geodesic_transport(t,Yold_tag,Ynew_tag)
258      implicit none
259      real*8   t
260      integer  Yold_tag
261      integer  Ynew_tag
262
263#include "bafdecls.fh"
264#include "errquit.fh"
265#include "c_geodesic_common.fh"
266
267*     **** local variables ****
268      complex*16 zero,one,mone
269      parameter (zero=(0.0d0,0.0d0), one=(1.0d0,0.0d0))
270      parameter (mone=(-1.0d0,0.0d0))
271
272      logical    value
273      integer    nb,npack1,nemax,nbrillq
274      integer    tmp1(2),tmp2(2),tmp3(2),tmpC(2),tmpS(2)
275      integer    yoldshift,ynewshift,ushift,sshift,vtshift
276
277*     **** external functions ****
278      logical  Pneb_w_push_get,Pneb_w_pop_stack
279      integer  cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
280      external Pneb_w_push_get,Pneb_w_pop_stack
281      external cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
282
283      call nwpw_timing_start(10)
284      call Cram_max_npack(npack1)
285      nemax   = cpsi_ne(1)+cpsi_ne(2)
286      nbrillq = cpsi_nbrillq()
287
288*     **** allocate tmp space ****
289      value =           Pneb_w_push_get(0,1,tmp1)
290      value = value.and.Pneb_w_push_get(0,1,tmp2)
291      value = value.and.Pneb_w_push_get(0,1,tmp3)
292      value = value.and.BA_push_get(mt_dbl,nemax,'tmpC',tmpC(2),tmpC(1))
293      value = value.and.BA_push_get(mt_dbl,nemax,'tmpS',tmpS(2),tmpS(1))
294      if (.not. value)
295     > call errquit('c_geodesic_transport:out of stack memory',0,MA_ERR)
296
297      do nb=1,nbrillq
298         yoldshift = cpsi_data_get_chnk(Yold_tag,nb)
299         ynewshift = cpsi_data_get_chnk(Ynew_tag,nb)
300         ushift    = cpsi_data_get_chnk(   U_tag,nb)
301         sshift    = cpsi_data_get_chnk(   S_tag,nb)
302         vtshift   = cpsi_data_get_chnk(  Vt_tag,nb)
303
304         call Pneb_SCVtrans2(0,nb,t,
305     >                   dbl_mb(sshift),
306     >                   dbl_mb(vtshift),
307     >                   dcpl_mb(tmp1(1)),
308     >                   dcpl_mb(tmp3(1)),
309     >                   dbl_mb(tmpC(1)),
310     >                   dbl_mb(tmpS(1)))
311         call Pneb_www_Multiply2(0,nb,
312     >                   one,
313     >                   dbl_mb(vtshift),
314     >                   dcpl_mb(tmp1(1)),
315     >                   zero,
316     >                   dcpl_mb(tmp2(1)))
317
318         call Pneb_fwf_Multiply(0,nb,
319     >                   mone,
320     >                   dbl_mb(yoldshift),npack1,
321     >                   dcpl_mb(tmp2(1)),
322     >                   zero,
323     >                   dbl_mb(ynewshift))
324
325         call Pneb_fwf_Multiply(0,nb,
326     >                   one,
327     >                   dbl_mb(ushift),npack1,
328     >                   dcpl_mb(tmp3(1)),
329     >                   one,
330     >                   dbl_mb(ynewshift))
331      end do
332*     **** deallocate tmp space ****
333      value =           BA_pop_stack(tmpS(2))
334      value = value.and.BA_pop_stack(tmpC(2))
335      value = value.and.Pneb_w_pop_stack(tmp3)
336      value = value.and.Pneb_w_pop_stack(tmp2)
337      value = value.and.Pneb_w_pop_stack(tmp1)
338      if (.not. value)
339     > call errquit('c_geodesic_transport:error popping stack',0,MA_ERR)
340
341      call nwpw_timing_end(10)
342
343      return
344      end
345
346
347*     ***********************************
348*     *					*
349*     *		c_geodesic_Gtransport	*
350*     *					*
351*     ***********************************
352*
353*     Uses - geodesic common block
354*
355
356      subroutine c_geodesic_Gtransport(t,Yold_tag,tG_tag)
357      implicit none
358      real*8   t
359      integer  Yold_tag
360      integer  tG_tag
361
362#include "bafdecls.fh"
363#include "errquit.fh"
364#include "c_geodesic_common.fh"
365
366*     **** local variables ****
367      complex*16 zero,one,mone
368      parameter (zero=(0.0d0,0.0d0), one=(1.0d0,0.0d0))
369      parameter (mone=(-1.0d0,0.0d0))
370
371      logical  value
372      integer  nb,npack1,nemax,nbrillq
373      integer  tmp1(2),tmp2(2),tmp3(2),tmpC(2),tmpS(2)
374      integer  yoldshift,tGshift,ushift,sshift,vtshift
375
376*     **** external functions ****
377      logical  Pneb_w_push_get,Pneb_w_pop_stack
378      integer  cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
379      external Pneb_w_push_get,Pneb_w_pop_stack
380      external cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
381
382      call nwpw_timing_start(10)
383      call Cram_max_npack(npack1)
384      nemax   = cpsi_ne(1)+cpsi_ne(2)
385      nbrillq = cpsi_nbrillq()
386
387
388*     **** allocate tmp space ****
389      value =           Pneb_w_push_get(0,1,tmp1)
390      value = value.and.Pneb_w_push_get(0,1,tmp2)
391      value = value.and.Pneb_w_push_get(0,1,tmp3)
392      value = value.and.BA_push_get(mt_dbl,nemax,'tmpC',tmpC(2),tmpC(1))
393      value = value.and.BA_push_get(mt_dbl,nemax,'tmpS',tmpS(2),tmpS(1))
394      if (.not. value)
395     > call errquit('c_geodesic_transport:out of stack memory',0,MA_ERR)
396
397
398      do nb=1,nbrillq
399         yoldshift = cpsi_data_get_chnk(Yold_tag,nb)
400         tGshift   = cpsi_data_get_chnk(  tG_tag,nb)
401         ushift    = cpsi_data_get_chnk(   U_tag,nb)
402         sshift    = cpsi_data_get_chnk(   S_tag,nb)
403         vtshift   = cpsi_data_get_chnk(  Vt_tag,nb)
404
405         call Pneb_ffw_Multiply(0,nb,
406     >                   dbl_mb(ushift),
407     >                   dbl_mb(tGshift),npack1,
408     >                   dcpl_mb(tmp2(1)))
409         call Pneb_SCVtrans3(0,nb,t,
410     >                   dbl_mb(sshift),
411     >                   dcpl_mb(tmp2(1)),
412     >                   dcpl_mb(tmp1(1)),
413     >                   dcpl_mb(tmp3(1)),
414     >                   dbl_mb(tmpC(1)),
415     >                   dbl_mb(tmpS(1)))
416         call Pneb_www_Multiply2(0,nb,
417     >                   one,
418     >                   dbl_mb(vtshift),
419     >                   dcpl_mb(tmp1(1)),
420     >                   zero,
421     >                   dcpl_mb(tmp2(1)))
422         call Pneb_fwf_Multiply(0,nb,
423     >                   mone,
424     >                   dbl_mb(yoldshift),npack1,
425     >                   dcpl_mb(tmp2(1)),
426     >                   one,
427     >                   dbl_mb(tGshift))
428         call Pneb_fwf_Multiply(0,nb,
429     >                   mone,
430     >                   dbl_mb(ushift),npack1,
431     >                   dcpl_mb(tmp3(1)),
432     >                   one,
433     >                   dbl_mb(tGshift))
434      end do
435
436*     **** deallocate tmp space ****
437      value =           BA_pop_stack(tmpS(2))
438      value = value.and.BA_pop_stack(tmpC(2))
439      value = value.and.Pneb_w_pop_stack(tmp3)
440      value = value.and.Pneb_w_pop_stack(tmp2)
441      value = value.and.Pneb_w_pop_stack(tmp1)
442      if (.not. value)
443     > call errquit('c_geodesic_transport:error popping stack',0,MA_ERR)
444
445      call nwpw_timing_end(10)
446
447      return
448      end
449
450
451*     *******************************************
452*     *						*
453*     *		c_geodesic_transport_junk	*
454*     *						*
455*     *******************************************
456*
457*   Temporary code until BGrsm_list fixed
458*     Uses - geodesic common block
459*
460
461      subroutine c_geodesic_transport_junk(t,Yold_tag,Ynew)
462      implicit none
463      real*8   t
464      integer Yold_tag
465      complex*16  Ynew(*)
466
467#include "bafdecls.fh"
468#include "errquit.fh"
469#include "c_geodesic_common.fh"
470
471*     **** local variables ****
472      complex*16 zero,one,mone
473      parameter (zero=(0.0d0,0.0d0), one=(1.0d0,0.0d0))
474      parameter (mone=(-1.0d0,0.0d0))
475
476      logical    value
477      integer    nb,npack1,nemax,nbrillq
478      integer    tmp1(2),tmp2(2),tmp3(2),tmpC(2),tmpS(2)
479c      integer    yoldshift,ynewshift,ushift,sshift,vtshift
480      integer    yoldshift,ushift,sshift,vtshift
481      integer    nbshift
482
483*     **** external functions ****
484      logical  Pneb_w_push_get,Pneb_w_pop_stack
485      integer  cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
486      external Pneb_w_push_get,Pneb_w_pop_stack
487      external cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
488
489      call nwpw_timing_start(10)
490      call Cram_max_npack(npack1)
491      nemax   = cpsi_ne(1)+cpsi_ne(2)
492      nbrillq = cpsi_nbrillq()
493      nbshift = nemax*npack1
494
495*     **** allocate tmp space ****
496      value =           Pneb_w_push_get(0,1,tmp1)
497      value = value.and.Pneb_w_push_get(0,1,tmp2)
498      value = value.and.Pneb_w_push_get(0,1,tmp3)
499      value = value.and.BA_push_get(mt_dbl,nemax,'tmpC',tmpC(2),tmpC(1))
500      value = value.and.BA_push_get(mt_dbl,nemax,'tmpS',tmpS(2),tmpS(1))
501      if (.not. value)
502     > call errquit('c_geodesic_transport:out of stack memory',0,MA_ERR)
503
504      do nb=1,nbrillq
505         yoldshift = cpsi_data_get_chnk(Yold_tag,nb)
506c         ynewshift = cpsi_data_get_chnk(Ynew_tag,nb)
507         ushift    = cpsi_data_get_chnk(   U_tag,nb)
508         sshift    = cpsi_data_get_chnk(   S_tag,nb)
509         vtshift   = cpsi_data_get_chnk(  Vt_tag,nb)
510
511         call Pneb_SCVtrans2(0,nb,t,
512     >                   dbl_mb(sshift),
513     >                   dbl_mb(vtshift),
514     >                   dcpl_mb(tmp1(1)),
515     >                   dcpl_mb(tmp3(1)),
516     >                   dbl_mb(tmpC(1)),
517     >                   dbl_mb(tmpS(1)))
518         call Pneb_www_Multiply2(0,nb,
519     >                   one,
520     >                   dbl_mb(vtshift),
521     >                   dcpl_mb(tmp1(1)),
522     >                   zero,
523     >                   dcpl_mb(tmp2(1)))
524
525         call Pneb_fwf_Multiply(0,nb,
526     >                   mone,
527     >                   dbl_mb(yoldshift),npack1,
528     >                   dcpl_mb(tmp2(1)),
529     >                   zero,
530     >                   Ynew(1+(nb-1)*nbshift))
531
532         call Pneb_fwf_Multiply(0,nb,
533     >                   one,
534     >                   dbl_mb(ushift),npack1,
535     >                   dcpl_mb(tmp3(1)),
536     >                   one,
537     >                   Ynew(1+(nb-1)*nbshift))
538      end do
539*     **** deallocate tmp space ****
540      value =           BA_pop_stack(tmpS(2))
541      value = value.and.BA_pop_stack(tmpC(2))
542      value = value.and.Pneb_w_pop_stack(tmp3)
543      value = value.and.Pneb_w_pop_stack(tmp2)
544      value = value.and.Pneb_w_pop_stack(tmp1)
545      if (.not. value)
546     > call errquit('c_geodesic_transport:error popping stack',0,MA_ERR)
547
548      call nwpw_timing_end(10)
549
550      return
551      end
552
553
554*     *******************************************
555*     *						*
556*     *		c_geodesic_Gtransport_junk	*
557*     *						*
558*     *******************************************
559*
560*   Temporary code until BGrsm_list fixed
561
562*     Uses - geodesic common block
563*
564
565      subroutine c_geodesic_Gtransport_junk(t,Yold_tag,tG)
566      implicit none
567      real*8   t
568      integer   Yold_tag
569      complex*16 tG(*)
570c      integer  tG_tag
571
572#include "bafdecls.fh"
573#include "errquit.fh"
574#include "c_geodesic_common.fh"
575
576*     **** local variables ****
577      complex*16 zero,one,mone
578      parameter (zero=(0.0d0,0.0d0), one=(1.0d0,0.0d0))
579      parameter (mone=(-1.0d0,0.0d0))
580
581      logical  value
582      integer  nb,npack1,nemax,nbrillq
583      integer  tmp1(2),tmp2(2),tmp3(2),tmpC(2),tmpS(2)
584c      integer  yoldshift,tGshift,ushift,sshift,vtshift
585      integer  yoldshift,ushift,sshift,vtshift
586      integer  nbshift
587
588*     **** external functions ****
589      logical  Pneb_w_push_get,Pneb_w_pop_stack
590      integer  cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
591      external Pneb_w_push_get,Pneb_w_pop_stack
592      external cpsi_ne,cpsi_nbrillq,cpsi_data_get_chnk
593
594      call nwpw_timing_start(10)
595      call Cram_max_npack(npack1)
596      nemax   = cpsi_ne(1)+cpsi_ne(2)
597      nbrillq = cpsi_nbrillq()
598      nbshift = nemax*npack1
599
600
601*     **** allocate tmp space ****
602      value =           Pneb_w_push_get(0,1,tmp1)
603      value = value.and.Pneb_w_push_get(0,1,tmp2)
604      value = value.and.Pneb_w_push_get(0,1,tmp3)
605      value = value.and.BA_push_get(mt_dbl,nemax,'tmpC',tmpC(2),tmpC(1))
606      value = value.and.BA_push_get(mt_dbl,nemax,'tmpS',tmpS(2),tmpS(1))
607      if (.not. value)
608     > call errquit('c_geodesic_transport:out of stack memory',0,MA_ERR)
609
610
611      do nb=1,nbrillq
612         yoldshift = cpsi_data_get_chnk(Yold_tag,nb)
613c         tGshift   = cpsi_data_get_chnk(  tG_tag,nb)
614         ushift    = cpsi_data_get_chnk(   U_tag,nb)
615         sshift    = cpsi_data_get_chnk(   S_tag,nb)
616         vtshift   = cpsi_data_get_chnk(  Vt_tag,nb)
617
618         call Pneb_ffw_Multiply(0,nb,
619     >                   dbl_mb(ushift),
620     >                   tG(1+(nb-1)*nbshift),npack1,
621     >                   dcpl_mb(tmp2(1)))
622         call Pneb_SCVtrans3(0,nb,t,
623     >                   dbl_mb(sshift),
624     >                   dcpl_mb(tmp2(1)),
625     >                   dcpl_mb(tmp1(1)),
626     >                   dcpl_mb(tmp3(1)),
627     >                   dbl_mb(tmpC(1)),
628     >                   dbl_mb(tmpS(1)))
629         call Pneb_www_Multiply2(0,nb,
630     >                   one,
631     >                   dbl_mb(vtshift),
632     >                   dcpl_mb(tmp1(1)),
633     >                   zero,
634     >                   dcpl_mb(tmp2(1)))
635         call Pneb_fwf_Multiply(0,nb,
636     >                   mone,
637     >                   dbl_mb(yoldshift),npack1,
638     >                   dcpl_mb(tmp2(1)),
639     >                   one,
640     >                   tG(1+(nb-1)*nbshift))
641         call Pneb_fwf_Multiply(0,nb,
642     >                   mone,
643     >                   dbl_mb(ushift),npack1,
644     >                   dcpl_mb(tmp3(1)),
645     >                   one,
646     >                   tG(1+(nb-1)*nbshift))
647      end do
648
649*     **** deallocate tmp space ****
650      value =           BA_pop_stack(tmpS(2))
651      value = value.and.BA_pop_stack(tmpC(2))
652      value = value.and.Pneb_w_pop_stack(tmp3)
653      value = value.and.Pneb_w_pop_stack(tmp2)
654      value = value.and.Pneb_w_pop_stack(tmp1)
655      if (.not. value)
656     > call errquit('c_geodesic_transport:error popping stack',0,MA_ERR)
657
658      call nwpw_timing_end(10)
659
660      return
661      end
662
663