ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
_sol.f
Go to the documentation of this file.
1  subroutine f_solve(isol,wdm)
2 
3 !--ISOL-->if ISOL equal 1 reodering and factorization are asumed
4 !-- to be done previosly.(PATH=3 will be used in SDRVD)
5 !--
6 !-- if ISOL equal 0 reodering and factorization will be done
7 !-- (PATH=1 will be used in ODRVD and SDRVD)
8 !--
9 !--WDM two dimensional grid array containing rezult of solving
10 
11  include 'double.inc'
12  include 'parevo.inc'
13  parameter(nkp=njlim)
14  include 'dim.inc'
15  include 'compol.inc'
16  include 'compol_add.inc'
17 
18  common /comwrc/ rsp,p,ip
19 
20  real*8 zw(neqp),rsp(nspp),wdm(nrp,ntp)
21 
22  integer p(neqp),ip(neqp),isp(nspp),ipath,flag,esp
23 
24  equivalence(rsp(1),isp(1))
25 
26 
27  ! write(*,*) 'solve:enter'
28 
29 
30  if( itin/nitdel*nitdel+nitbeg .eq. itin
31  * .OR. itin.lt.nitbeg ) then
32 
33 
34  ipath=3
35 
36  if(isol.ne.0) go to 20
37 
38  call odrvd(neq,ia,ja,a,p,ip,nspp,isp,1,flag)
39 
40  !do 10 i=1,neqp
41  ! ip(i)=i
42  ! p(i)=i
43 !10 continue
44 
45  !write(6,*) 'odrv flag=',flag
46 
47  ipath=1
48  !isol=1
49 
50  20 continue
51 
52 
53  call sdrvd(neq,p,ip,ia,ja,a,right,zw,nspp,
54  * isp,rsp,esp,ipath,flag)
55 
56 c do 860 i=1,neq
57 c write(6,*) 'zw(i) i',i,zw(i)
58 c860 continue
59 
60  ! write(6,*) 'sdrv: flag,esp',flag,esp
61 
62  !call f_nev(zw)
63 
64 
65  else
66 
67  call solvit(isol,zw)
68 
69  endif
70 
71 c raspakovka reshenia
72 
73  do 500 i=1,nr1
74  do 500 j=2,nt1
75 
76  ieq=numlin(i,j,nr,nt)
77 
78  wdm(i,j)=zw(ieq)
79 
80  500 continue
81 
82  do 600 i=1,nr
83 
84  wdm(i,1)=wdm(i,nt1)
85  wdm(i,nt)=wdm(i,2)
86 
87  600 continue
88 
89  return
90  end
91 
92 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
93 
94  subroutine solvit(isol,zw)
95 
96  include 'double.inc'
97  include 'parevo.inc'
98  parameter(nkp=njlim)
99  include 'dim.inc'
100  include 'compol.inc'
101  include 'compol_add.inc'
102 
103  common /comwrc/ rsp,p,ip
104 
105  real*8 zw(neqp),zyy(neqp),rsp(nspp),wpp(neqp),wzz(neqp),
106  * wrr(neqp),zuu(neqp)
107 
108  integer p(neqp),ip(neqp),isp(nspp),ipath,flag,esp
109 
110  equivalence(rsp(1),isp(1))
111 
112  !write(6,*) 'solvit:enter'
113 
114  itmax=10
115  relerr=1.d-8
116 
117  do il=1,neq
118  wrr(il)=right(il)
119 
120  enddo
121 
122  if(isol.eq.0) then
123 
124  do i=1,nr1
125  do j=2,nt1
126 
127  ieq=numlin(i,j,nr,nt)
128 
129  zw(ieq)=g(i,j)
130 
131  enddo
132  enddo
133 
134  elseif(isol.eq.1) then
135 
136  do i=1,nr1
137  do j=2,nt1
138 
139  ieq=numlin(i,j,nr,nt)
140 
141  zw(ieq)=psii(i,j)
142 
143  enddo
144  enddo
145  endif
146 
147  do il=1,neq
148 
149  i1=ia(il)
150  i2=ia(il+1)-1
151 
152  znes=0.d0
153 
154  do im=i1,i2
155 
156  ic=ja(im)
157 
158  znes=znes+daop(im)*zw(ic)
159 
160  enddo
161 
162  zyy(il)=wrr(il)-znes
163 
164  enddo
165 
166  call sdrvd(neq,p,ip,ia,ja,aop0,zyy,zw,nspp,
167  * isp,rsp,esp,3,flag)
168 
169  !if(itin/2*2.ne.itin)
170  return
171 
172  itk=0
173  ido=0
174 
175  77 call dpcgrc(ido,neq,zw,wpp,wrr,wzz,relerr,itmax)
176 
177  !write(6,*) 'ido',ido
178  ! pause ' '
179 
180  if(ido.eq.1) then
181 
182  do 10 il=1,neq
183 
184  i1=ia(il)
185  i2=ia(il+1)-1
186 
187  znes=0.d0
188 
189  do 100 im=i1,i2
190 
191  ic=ja(im)
192 
193  znes=znes+a(im)*wpp(ic)
194 
195  100 continue
196 
197  wzz(il)=znes
198 
199  10 continue
200 
201  go to 77
202 
203  elseif(ido.eq.2) then
204 
205  call sdrvd(neq,p,ip,ia,ja,aop0,wrr,wzz,nspp,
206  * isp,rsp,esp,3,flag)
207  ! write(6,*) 'sdrv: flag,esp',flag,esp
208 
209  if(itk.eq.1) go to 1054
210 
211  do il=1,neq
212 
213  i1=ia(il)
214  i2=ia(il+1)-1
215 
216  znes=0.d0
217 
218  do im=i1,i2
219 
220  ic=ja(im)
221 
222  znes=znes+daop(im)*wzz(ic)
223 
224  enddo
225 
226  zyy(il)=znes
227 
228  enddo
229 
230  call sdrvd(neq,p,ip,ia,ja,aop0,zyy,zuu,nspp,
231  * isp,rsp,esp,3,flag)
232  !write(6,*) 'sdrv: flag,esp',flag,esp
233 
234  do il=1,neq
235 
236  wzz(il)=wzz(il)-zuu(il)
237 
238  enddo
239 
240  1054 continue
241 
242  ! zermax=0.d0
243  ! znvmax=0.d0
244  ! do il=1,neq
245 
246  ! znev=dabs(wrr(il))
247  ! znvmax=dmax1(znev,znvmax)
248 
249  !zerr=dabs(zw(il)-zw0(il))
250  !zermax=dmax1(zerr,zermax)
251 
252  ! enddo
253 
254  itk=itk + 1
255  ! write(6,*) 'nev,err iter ',znvmax,zermax,itk
256 
257 
258  ! if(znvmax.lt.relerr) go to 537
259  go to 77
260 
261  endif
262 
263  537 continue
264 
265  !write(6,*) 'nev,err iter ',znvmax,zermax,itk
266  !! write(6,*) ' itk ',itk
267 
268  ! pause ' '
269  ! stop
270  return
271  end
272 
273 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
274 
275  subroutine f_nev(zw)
276 
277  include 'double.inc'
278  include 'parevo.inc'
279  parameter(nkp=njlim)
280  include 'dim.inc'
281  include 'compol.inc'
282  include 'compol_add.inc'
283 
284  real*8 zw(neqp)
285 
286  do 10 il=1,neq
287 
288  i1=ia(il)
289  i2=ia(il+1)-1
290 
291  znev=0.d0
292  znmx=0.d0
293 
294  do 100 im=i1,i2
295 
296  ic=ja(im)
297 
298  znev=znev+a(im)*zw(ic)
299 
300  ! if(il.eq.1) then
301  !
302  ! write(6,*) 'im a(im)',im,a(im)
303  ! write(6,*) 'ic zw(ic)',ic,zw(ic)
304  ! write(6,*) 'znev ' ,znev
305  ! endif
306 
307  100 continue
308 
309  znev=right(il)-znev
310  znab=dabs(znev)
311 c write(6,*) '***right znev ***',il,right(il),znev
312  znmx=dmax1(znmx,znab)
313 
314  10 continue
315  !write(6,*) 'nev:',znmx
316 
317  return
318  end
319 
320 
321 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
322 
323  subroutine f_solint
324 
325  include 'double.inc'
326  include 'parevo.inc'
327  parameter(nkp=njlim)
328  include 'dim.inc'
329  include 'compol.inc'
330  include 'compol_add.inc'
331 
332  common /comwrp/ rsp1,p1,ip1
333 
334  real*8 zw(neqp),rsp1(nspp)
335 
336  integer p1(neqp),ip1(neqp),isp1(nspp),ipath,flag,esp
337 
338  equivalence(rsp1(1),isp1(1))
339  common
340  * /c_kpr/kpr
341 
342  ! write(6,*) 'solint:enter'
343 
344  if( itin/nitdel*nitdel+nitbeg .eq. itin
345  * .OR. itin.lt.nitbeg ) then
346 
347  call odrvd(neqpla,ia,ja,a,p1,ip1,nspp,isp1,1,flag)
348 
349  !do 10 i=1,neqpla
350  ! ip1(i)=i
351  ! p1(i)=i
352 ! 10 continue
353 
354  ! write(6,*) 'odrv flag=',flag
355 
356  ipath=1
357 
358  20 continue
359 
360  call sdrvd(neqpla,p1,ip1,ia,ja,a,right,zw,nspp,
361  * isp1,rsp1,esp,ipath,flag)
362 
363 c do 860 i=1,neq
364 c write(6,*) 'zw(i) i',i,zw(i)
365 c860 continue
366 
367  ! write(6,*) 'sdrv: flag,esp',flag,esp
368  else
369 
370  call soleit(zw)
371 
372  endif
373 
374 c raspakovka reshenia
375  errpss=0.d0
376 
377  do 500 i=1,iplas-1
378  do 500 j=2,nt1
379 
380  ieq=numlin(i,j,nr,nt)
381 
382  delpsi=dabs(psi(i,j)-zw(ieq))
383  errpss=dmax1(errpss,delpsi)
384 
385  psi(i,j)=zw(ieq)
386 
387  500 continue
388 
389  do 600 i=1,iplas1
390 
391  psi(i,1)=psi(i,nt1)
392  psi(i,nt)=psi(i,2)
393 
394  600 continue
395 
396 
397  if(kpr.eq.1) then
398  write(*,*) 'solint:errpss',errpss
399  endif
400 
401 
402  return
403  end
404 
405 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
406 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
407 
408 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
409 
410  subroutine solext
411 
412  include 'double.inc'
413  include 'parevo.inc'
414  parameter(nkp=njlim)
415  include 'dim.inc'
416  include 'compol.inc'
417  include 'compol_add.inc'
418 
419  common /comwrp/ rsp1,p1,ip1
420 
421  real*8 zw(neqp),rsp1(nspp)
422 
423  integer p1(neqp),ip1(neqp),isp1(nspp),ipath,flag,esp
424 
425  equivalence(rsp1(1),isp1(1))
426 
427  write(*,*) 'solext:enter'
428 
429  if( itin/nitdel*nitdel+nitbeg .eq. itin
430  * .OR. itin.lt.nitbeg ) then
431 
432  call odrvd(neqpla,ia,ja,a,p1,ip1,nspp,isp1,1,flag)
433 
434  !do 10 i=1,neqpla
435  ! ip1(i)=i
436  ! p1(i)=i
437 ! 10 continue
438 
439  !write(6,*) 'odrv flag=',flag
440 
441  ipath=1
442 
443  20 continue
444 
445  call sdrvd(neqpla,p1,ip1,ia,ja,a,right,zw,nspp,
446  * isp1,rsp1,esp,ipath,flag)
447 
448 c do 860 i=1,neq
449 c write(6,*) 'zw(i) i',i,zw(i)
450 c860 continue
451 
452  write(*,*) 'sdrv: flag,esp',flag,esp
453  else
454 
455  call soleit(zw)
456 
457  endif
458 
459 c raspakovka reshenia
460 
461  do 500 i=1,iplas-1
462  do 500 j=2,nt1
463 
464  ieq=numlin(i,j,nr,nt)
465 
466  psie(i,j)=zw(ieq)
467 
468  500 continue
469 
470  do 600 i=1,nr
471 
472  psie(i,1)=psie(i,nt1)
473  psie(i,nt)=psie(i,2)
474 
475  600 continue
476 
477  return
478  end
479 
480 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
482 
483  subroutine soleit(zw)
484 
485  include 'double.inc'
486  include 'parevo.inc'
487  parameter(nkp=njlim)
488  include 'dim.inc'
489  include 'compol.inc'
490  include 'compol_add.inc'
491 
492  common /comwrp/ rsp1,p1,ip1
493 
494  real*8 zw(neqp),zyy(neqp),rsp1(nspp),wpp(neqp),wzz(neqp),
495  * wrr(neqp),zuu(neqp)
496 
497  integer p1(neqp),ip1(neqp),isp1(nspp),ipath,flag,esp
498 
499  equivalence(rsp1(1),isp1(1))
500 
501  !write(6,*) 'soleit:enter'
502 
503  itmax=10
504  relerr=1.d-8
505 
506 
507  do i=1,nr1
508  do j=2,nt1
509 
510  ieq=numlin(i,j,nr,nt)
511 
512  zw(ieq)=psie(i,j)
513 
514  enddo
515  enddo
516 
517  do il=1,neqpla
518  wrr(il)=right(il)
519 
520  enddo
521 
522  do il=1,neqpla
523 
524  i1=ia(il)
525  i2=ia(il+1)-1
526 
527  znes=0.d0
528 
529  do im=i1,i2
530 
531  ic=ja(im)
532 
533  znes=znes+dapp(im)*zw(ic)
534 
535  enddo
536 
537  zyy(il)=wrr(il)-znes
538 
539  enddo
540 
541  call sdrvd(neqpla,p1,ip1,ia,ja,app0,zyy,zw,nspp,
542  * isp1,rsp1,esp,3,flag)
543 
544  !if(itin/2*2.ne.itin)
545  return
546 
547  itk=0
548  ido=0
549 
550  77 call dpcgrc(ido,neqpla,zw,wpp,wrr,wzz,relerr,itmax)
551 
552  !write(6,*) 'ido',ido
553  ! pause 'soleit '
554 
555  if(ido.eq.1) then
556 
557  do 10 il=1,neqpla
558 
559  i1=ia(il)
560  i2=ia(il+1)-1
561 
562  znes=0.d0
563 
564  do 100 im=i1,i2
565 
566  ic=ja(im)
567 
568  znes=znes+a(im)*wpp(ic)
569 
570  100 continue
571 
572  wzz(il)=znes
573 
574  10 continue
575 
576  go to 77
577 
578  elseif(ido.eq.2) then
579 
580  call sdrvd(neqpla,p1,ip1,ia,ja,app0,wrr,wzz,nspp,
581  * isp1,rsp1,esp,3,flag)
582  ! write(6,*) 'sdrv: flag,esp',flag,esp
583 
584  if(itk.eq.1) go to 1056
585 
586  do il=1,neqpla
587 
588  i1=ia(il)
589  i2=ia(il+1)-1
590 
591  znes=0.d0
592 
593  do im=i1,i2
594 
595  ic=ja(im)
596 
597  znes=znes+dapp(im)*wzz(ic)
598 
599  enddo
600 
601  zyy(il)=znes
602 
603  enddo
604 
605  call sdrvd(neqpla,p1,ip1,ia,ja,app0,zyy,zuu,nspp,
606  * isp1,rsp1,esp,3,flag)
607  ! write(6,*) 'sdrv: flag,esp',flag,esp
608 
609  do il=1,neqpla
610 
611  wzz(il)=wzz(il)-zuu(il)
612 
613  enddo
614 
615  1056 continue
616 
617 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
618 ! do il=1,neq
619 !
620 ! i1=ia(il)
621 ! i2=ia(il+1)-1
622 !
623 ! znes=0.d0
624 !
625 ! do im=i1,i2
626 !
627 ! ic=ja(im)
628 !
629 ! if(ic.eq.il) then
630 ! adiag=a(im)
631 ! go to 110
632 ! endif
633 !
634 ! enddo
635 !
636 ! 110 wzz(il)=wrr(il)/adiag
637 !
638 ! enddo
639 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
640 
641  ! zermax=0.d0
642  ! znvmax=0.d0
643  ! do il=1,neq
644  !
645  ! znev=dabs(wrr(il))
646  ! znvmax=dmax1(znev,znvmax)
647 
648  ! zerr=dabs(zw(il)-zw0(il))
649  ! zermax=dmax1(zerr,zermax)
650 
651  ! enddo
652 
653  itk=itk + 1
654  ! write(6,*) 'nev,err iter ',znvmax,zermax,itk
655 
656 
657  ! if(znvmax.lt.relerr) go to 537
658  go to 77
659 
660  endif
661 
662  537 continue
663 
664  !write(6,*) 'nev,err iter ',znvmax,zermax,itk
665  !write(6,*) ' iter ',itk
666 
667  ! pause ' '
668  ! stop
669  return
670  end
671 
672 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
673 
674 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
675  subroutine ctra
676 
677  include'double.inc'
678  include 'parevo.inc'
679  parameter(nkp=njlim)
680  include'dim.inc'
681  include'compol.inc'
682  include 'compol_add.inc'
683 
684 c real*8 a(lp)
685  real*8 zw(neqp),ssw(neqp)
686 c equivalence(a(1),vol1(1,1))
687 
688  sqcen=0.d0
689  do 30 j=2,nt1
690 
691  sqcen=sqcen+sq1(1,j)+sq4(1,j)
692 
693  30 continue
694  ssw(1)=sqcen
695 
696  do 20 i=2,nr1
697  do 20 j=2,nt1
698 
699  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
700 
701  il=numlin(i,j,nr,nt)
702  ssw(il)=sqk
703  20 continue
704 
705  do 500 i=1,nr1
706  do 500 j=2,nt1
707 
708  il=numlin(i,j,nr,nt)
709 
710  !zw(il)=1.
711  !zw(il)=ro(i,j)**2
712  !zw(il)=z(i,j)
713  !zw(il)=r(i,j)**2
714  zw(il)=psie(i,j)
715 
716  500 continue
717 
718 
719 
720  znmx=0.d0
721  do 10 il=1,neq
722 
723  i1=ia(il)
724  i2=ia(il+1)-1
725 
726  znev=0.d0
727 
728  do 100 im=i1,i2
729 
730  ic=ja(im)
731 
732  znev=znev+a(im)*zw(ic)
733 
734  ! if(il.eq.1) then
735  !
736  ! write(6,*) 'im a(im)',im,a(im)
737  ! write(6,*) 'ic zw(ic)',ic,zw(ic)
738  ! write(6,*) 'znev ' ,znev
739  ! endif
740 
741  100 continue
742  znev=znev/ssw(il)
743  znab=dabs(znev)
744  !write(1,*) '*** znev ***',il,znev
745  znmx=dmax1(znmx,znab)
746 
747  10 continue
748 
749  !write(1,*) 'nev:',znmx
750 
751  return
752  end
753 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
754  subroutine aprod
755 
756  include 'double.inc'
757  include 'parevo.inc'
758  parameter(nkp=njlim)
759  include 'dim.inc'
760  include 'compol.inc'
761  include 'compol_add.inc'
762 
763  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
764  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
765 
766 c real*8 a(lp)
767  real*8 zw(nrp,ntp)
768 c equivalence(a(1),vol1(1,1))
769 
770  do 500 i=1,nr
771  do 500 j=1,nt
772 
773  zw(i,j)=psie(i,j)+clr*r(i,j)**2+clz*z(i,j)
774  !zw(i,j)= clr*r(i,j)**2+clz*z(i,j)
775  !zw(i,j)= r(i,j)**2
776  !zw(i,j)=psie(i,j)
777  !zw(i,j)=greeni(0.5d0,0.d0,r(i,j),z(i,j))
778 
779  500 continue
780 
781  do 10 i=1,nr1
782 
783  if(i.eq.1) then
784 
785  asum=0.d0
786  a0=0.d0
787 
788  do 20 j=2,nt1
789 
790  aj=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
791  a0=a0-aj
792  asum=asum+aj*zw(2,j)
793 
794  20 continue
795 
796  asum=asum+a0*zw(1,2)
797 
798  zpro(1)=asum
799 
800  else
801 
802  do 30 j=2,nt1
803 
804  a1=a13(i-1,j-1)
805  a2=a34(i-1,j-1)+a12(i-1,j)
806  a3=a24(i-1,j)
807  a4=a23(i-1,j-1)+a14(i,j-1)
808  a6=a14(i,j)+a23(i-1,j)
809  a7=a24(i,j-1)
810  a8=a34(i,j-1)+a12(i,j)
811  a9=a13(i,j)
812 
813  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
814 
815  u1=zw(i-1,j-1)
816  u2=zw(i-1,j)
817  u3=zw(i-1,j+1)
818  u4=zw(i,j-1)
819  u5=zw(i,j)
820  u6=zw(i,j+1)
821  u7=zw(i+1,j-1)
822  u8=zw(i+1,j)
823  u9=zw(i+1,j+1)
824 
825  il=numlin(i,j,nr,nt)
826 
827  zpro(il)=a1*u1+a2*u2+a3*u3+a4*u4+a5*u5+a6*u6+a7*u7+a8*u8+a9*u9
828 
829  30 continue
830 
831  endif
832 
833  10 continue
834 
835  do 200 i=1,nr1
836  do 200 j=2,nt1
837 
838  ieq=numlin(i,j,nr,nt)
839 
840  aex(i,j)=zpro(ieq)
841 
842  200 continue
843 
844  do 600 i=1,nr1
845 
846  aex(i,1)=aex(i,nt1)
847  aex(i,nt)=aex(i,2)
848 
849  600 continue
850 
851  return
852  end
853 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
854  subroutine aprod0
855 
856  include 'double.inc'
857  include 'parevo.inc'
858  parameter(nkp=njlim)
859  include 'dim.inc'
860  include 'compol.inc'
861  include 'compol_add.inc'
862 
863 c real*8 a(lp)
864  real*8 zw(neqp)
865 c equivalence(a(1),vol1(1,1))
866 
867 
868  do 500 i=1,nr1
869  do 500 j=2,nt1
870 
871  il=numlin(i,j,nr,nt)
872 
873  zw(il)=psie(i,j)+clr*r(i,j)**2+clz*z(i,j)
874  !zw(il)= clr*r(i,j)**2+clz*z(i,j)
875  !zw(il)= r(i,j)**2
876  !zw(il)=psie(i,j)
877  !zw(il)=greeni(0.5d0,0.d0,r(i,j),z(i,j))
878 
879  500 continue
880 
881  do 10 il=1,neq
882 
883  i1=ia(il)
884  i2=ia(il+1)-1
885 
886  asum=0.d0
887 
888  do 100 im=i1,i2
889  ic=ja(im)
890  asum=asum+a(im)*zw(ic)
891  100 continue
892 
893  zpro(il)=asum
894 
895  10 continue
896 
897  do 50 j=2,nt1
898  il=numlin(nr1,j,nr,nt)
899  zpro(il)=0.d0
900  50 continue
901 
902  do 200 i=1,nr1
903  do 200 j=2,nt1
904 
905  ieq=numlin(i,j,nr,nt)
906 
907  aex(i,j)=zpro(ieq)
908 
909  200 continue
910 
911  do 600 i=1,nr1
912 
913  aex(i,1)=aex(i,nt1)
914  aex(i,nt)=aex(i,2)
915 
916  600 continue
917 
918  return
919  end
920 
921 
aprod
subroutine aprod
Definition: _sol.f:754
odrvd
subroutine odrvd
Definition: SPAR.f:9
solext
subroutine solext
Definition: _sol.f:410
soleit
subroutine soleit(zw)
Definition: _sol.f:483
aprod0
subroutine aprod0
Definition: _sol.f:854
f_solint
subroutine f_solint
Definition: _sol.f:323
numlin
function numlin(i, j, nr, nt)
Definition: com_sub.f:1040
f_nev
subroutine f_nev(zw)
Definition: _sol.f:275
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
dpcgrc
subroutine dpcgrc(IDO, N, X, P, R, Z, RELERR, ITMAX)
Definition: DPCGRCV.f:1
ctra
subroutine ctra
Definition: _sol.f:675
sdrvd
subroutine sdrvd
Definition: SPAR.f:1411
solvit
subroutine solvit(isol, zw)
Definition: _sol.f:94
f_solve
subroutine f_solve(isol, wdm)
Definition: _sol.f:1