V3FIT
scalpot.f
1  SUBROUTINE scalpot(bvec, amatrix, ivacskip)
2  USE vacmod
3  USE parallel_include_module
4  USE timer_sub
5  IMPLICIT NONE
6 C-----------------------------------------------
7 C D u m m y A r g u m e n t s
8 C-----------------------------------------------
9  INTEGER, INTENT(in) :: ivacskip
10  REAL(dp), INTENT(out) :: bvec(mnpd2), amatrix(mnpd2*mnpd2)
11 C-----------------------------------------------
12 C L o c a l V a r i a b l e s
13 C-----------------------------------------------
14  INTEGER :: ip, istore, istart, istore_max, ndim
15  REAL(dp), ALLOCATABLE :: grpmn(:), green(:), gstore(:)
16  REAL(dp), ALLOCATABLE :: greenp(:,:)
17  REAL(dp) :: ton, toff, tonscal
18 C-----------------------------------------------
19  CALL second0(tonscal)
20 
21  IF (.NOT.ALLOCATED(amatsav)) THEN
22  stop 'AMATSAV: Allocation error in scalpot'
23  END IF
24 
25  ALLOCATE (grpmn(nuv3*mnpd2), stat=ip)
26  IF (ip .NE. 0) stop 'GRPMN: Allocation error in scalpot'
27 
28 !
29 ! COMPUTE TRANFORM OF ANALYTIC SOURCE AND KERNEL
30 ! ON EXIT, BVEC CONTAINS THE TRANSFORM OF THE ANALYTIC SOURCE
31 ! AND GRPMN CONTAINS TRANSFORM OF NORMAL DERIVATIVE
32 ! OF THE GREENS FUNCTION [PKM, EQ.(2.15)]
33 !
34 ! FOR ivacskip != 0, USE PREVIOUSLY COMPUTED bvecsav FOR SPEED
35 !
36 
37  ndim = mnpd2/mnpd
38  CALL analyt (grpmn, bvec, ivacskip, ndim)
39 
40  IF (ivacskip .NE. 0) THEN
41  bvec = bvec + bvecsav
42  ELSE
43 
44  istore_max = min(64,nuv3)
45 
46  ALLOCATE (green(nuv), gstore(nuv), greenp(nuv,istore_max),
47  & stat=ip)
48  IF (ip .NE. 0) stop 'Allocation error in scalpot'
49 ! bvecsav = bvec !Save in fouri now
50  gstore = 0
51 !
52 ! COMPUTE SURFACE INTEGRALS OF SOURCE AND GREENS FUNCTION NEEDED
53 ! FOR SPECTRAL DECOMPOSITION OF POTENTIAL INTEGRAL EQUATION
54 ! NOTE: SOURCE IS THE RHS OF EQ.(3.2), KERNEL IS THE LHS OF EQ (3.2).
55 ! IP IS THE INDEX OF THE PRIMED VARIABLE MESH.
56 !
57  istart = nuv3min - 1
58 
59  !SKS: Have to parallelize over ip since arrays computed in
60  !surface.f like rzb2, z1b, etc but used in green.f are
61  !known only within [nuv3min, nuv3max].
62 
63  primed: DO ip = nuv3min, nuv3max
64  istore = 1 + mod(ip-nuv3min,istore_max)
65 !
66 ! COMPUTE DIFFERENCE BETWEEN THE EXACT AND ANALYTIC GREENS FUNCTION AND GRADIENT
67 ! [FIRST TERMS IN EQ.(2.14, 2.16)].
68 !
69 ! BECAUSE OF THE LARGE SIZES INVOLVED (nuv3*nuv3), THIS IS DONE HERE BY STORING
70 ! THESE QUANTITIES - FOR ALL VALUES OF THE UNPRIMED U,V COORDINATE MESH - ON A
71 ! LIMITED SUBSET (ISTORE_max) OF PRIMED MESH VALUES (INDEXED BY ISTORE~IP),
72 ! THE FOURIER TRANSFORM OVER THE PRIMED MESH IS "BUILT-UP" BY MULTIPLE CALLS TO FOURP
73 ! WITHIN THIS LOOP.
74 !
75  CALL greenf (green, greenp(1,istore), ip)
76 
77 
78 ! PERFORM INTEGRAL (SUM) OVER PRIMED MESH OF NON-SINGULAR SOURCE TERM
79 ! [(h-hsing)(u,v,u',v') == bexni(ip)*green(u,v; ip) in Eq. 2.16]
80 ! AND STORE IT - FOR UNPRIMED MESH VALUES - IN GSTORE
81 
82  gstore = gstore + bexni(ip)*green
83 
84 !
85 ! PERFORM FOURIER INTEGRAL OF GRADIENT KERNEL (GREENP) OVER THE UNPRIMED MESH
86 ! AND STORE IN GRPMN (NOTE THAT GRPMN IS ADDED TO THE ANALYTIC PIECE IN EQ. 2.14,
87 ! - COMPUTED IN ANALYT - WHICH HAS THE APPROPRIATE SIN, COS FACTORS ALREADY)
88 !
89 
90  IF (istore.EQ.istore_max .OR. ip.EQ.nuv3max) THEN
91  CALL fourp (grpmn, greenp, istore, istart, ip, ndim)
92  END IF
93 
94  END DO primed
95 
96  CALL second0(ton)
97  IF (vlactive) THEN
98  CALL mpi_allreduce(mpi_in_place, gstore, SIZE(gstore),
99  & mpi_real8, mpi_sum, vac_comm, mpi_err)
100  END IF
101  CALL second0(toff)
102  allreduce_time = allreduce_time + (toff - ton)
103  timer_vac(tallr) = timer_vac(tallr) + (toff-ton)
104 !
105 ! COMPUTE FOURIER INTEGRAL OF GRADIENT (GRPMN) OVER PRIMED MESH IN EQ. 2.14
106 ! AND SOURCE (GSTORE) OVER UNPRIMED MESH IN EQ. 2.16
107 !
108  CALL fouri (grpmn, gstore, amatrix, amatsav, bvec,
109  & bvecsav, ndim)
110  DEALLOCATE (green, greenp, gstore)
111 
112  END IF
113 
114  DEALLOCATE (grpmn, stat=ip)
115 
116  amatrix = amatsav
117 
118 !
119 ! FINAL REDUCTION OVER VAC_COMM DONE ONCE HERE
120 !
121  CALL second0(ton)
122  IF (vlactive) THEN
123  CALL mpi_allreduce(mpi_in_place, bvec, SIZE(bvec), mpi_real8,
124  & mpi_sum, vac_comm, mpi_err)
125  END IF
126 
127  CALL second0(toff)
128  allreduce_time = allreduce_time + (toff - ton)
129  timer_vac(tanar) = timer_vac(tanar) + (toff-ton)
130 
131  scalpot_time = scalpot_time + (tonscal - toff)
132 
133  END SUBROUTINE scalpot