2 ! { dg-options "-O -fbounds-check -w" }
4 INTEGER, PARAMETER :: dp
= SELECTED_REAL_KIND ( 14, 200 )
5 INTEGER, DIMENSION(:), ALLOCATABLE
:: nco
,ncoset
,nso
,nsoset
6 INTEGER, DIMENSION(:,:,:), ALLOCATABLE
:: co
,coset
11 SUBROUTINE cossin(la_max
,npgfa
,zeta
,rpgfa
,la_min
,&
12 lb_max
,npgfb
,zetb
,rpgfb
,lb_min
,&
13 rac
,rbc
,kvec
,cosab
,sinab
)
14 REAL(KIND
=dp
), DIMENSION(ncoset(la_max
),&
15 ncoset(lb_max
)) :: sc
, ss
22 sc(coset(ax
,ay
,az
),1) = rap(1)*sc(coset(ax
-1,ay
,az
),1) +&
23 f2
* kvec(1)*ss(coset(ax
-1,ay
,az
),1)
24 ss(coset(ax
,ay
,az
),1) = rap(1)*ss(coset(ax
-1,ay
,az
),1) +&
25 f2
* kvec(1)*sc(coset(ax
-1,ay
,az
),1)
31 ss(1,coset(0,0,lb
)) = rbp(3)*ss(1,coset(0,0,lb
-1)) +&
32 f2
* kvec(3)*sc(1,coset(0,0,lb
-1))
35 ss(1,coset(bx
,by
,bz
)) = rbp(1)*ss(1,coset(bx
-1,by
,bz
)) +&
36 f2
* kvec(1)*sc(1,coset(bx
-1,by
,bz
))
42 DO j
=ncoset(lb_min
-1)+1,ncoset(lb_max
)
47 SUBROUTINE moment(la_max
,npgfa
,zeta
,rpgfa
,la_min
,&
48 lb_max
,npgfb
,zetb
,rpgfb
,&
50 REAL(KIND
=dp
), DIMENSION(:), INTENT(IN
) :: zeta
, rpgfa
51 REAL(KIND
=dp
), DIMENSION(:), INTENT(IN
) :: zetb
, rpgfb
52 REAL(KIND
=dp
), DIMENSION(:, :, :), &
54 REAL(KIND
=dp
), DIMENSION(3) :: rab
, rap
, rbp
, rpc
55 REAL(KIND
=dp
), DIMENSION(ncoset(la_max
),&
56 ncoset(lb_max
), ncoset(lc_max
)) :: s
59 IF (rpgfa(ipgf
) + rpgfb(jpgf
) < dab
) THEN
60 DO k
=1, ncoset(lc_max
)-1
61 DO j
=nb
+1,nb
+ncoset(lb_max
)
62 DO i
=na
+1,na
+ncoset(la_max
)
68 rpc
= zetp
*(zeta(ipgf
)*rac
+zetb(jpgf
)*rbc
)
69 DO l
=2, ncoset(lc_max
)
73 IF ( lz
> 1 ) l2
= coset(lx
,ly
,lz
-2)
74 ELSE IF ( ly
> 0 ) THEN
75 IF ( ly
> 1 ) l2
= coset(lx
,ly
-2,lz
)
76 IF ( lx
> 1 ) l2
= coset(lx
-2,ly
,lz
)
78 s(1,1,l
) = rpc(i
)*s(1,1,l1
)
79 IF ( l2
> 0 ) s(1,1,l
) = s(1,1,l
) + f2
*REAL(ni
,dp
)*s(1,1,l2
)
81 DO l
= 1, ncoset(lc_max
)
83 lx1
= coset(lx
-1,ly
,lz
)
86 ly1
= coset(lx
,ly
-1,lz
)
90 IF ( lz1
> 0 ) s(coset(0,0,la
),1,l
) = s(coset(0,0,la
),1,l
) + &
91 f2z
*s(coset(0,0,la
-1),1,lz1
)
92 IF ( ly1
> 0 ) s(coset(0,1,az
),1,l
) = s(coset(0,1,az
),1,l
) + &
93 f2y
*s(coset(0,0,az
),1,ly1
)
95 s(coset(0,ay
,az
),1,l
) = rap(2)*s(coset(0,ay
-1,az
),1,l
) +&
96 f2
*REAL(ay
-1,dp
)*s(coset(0,ay
-2,az
),1,l
)
97 IF ( ly1
> 0 ) s(coset(0,ay
,az
),1,l
) = s(coset(0,ay
,az
),1,l
) + &
98 f2y
*s(coset(0,ay
-1,az
),1,ly1
)
101 IF ( lx1
> 0 ) s(coset(1,ay
,az
),1,l
) = s(coset(1,ay
,az
),1,l
) + &
102 f2x
*s(coset(0,ay
,az
),1,lx1
)
106 s(coset(ax
,ay
,az
),1,l
) = rap(1)*s(coset(ax
-1,ay
,az
),1,l
) +&
107 f3
*s(coset(ax
-2,ay
,az
),1,l
)
108 IF ( lx1
> 0 ) s(coset(ax
,ay
,az
),1,l
) = s(coset(ax
,ay
,az
),1,l
) + &
109 f2x
*s(coset(ax
-1,ay
,az
),1,lx1
)
114 DO j
=2,ncoset(lb_max
)
115 DO i
=1,ncoset(la_max
)
119 DO la
=la_start
,la_max
-1
122 s(coset(ax
,ay
,az
),2,l
) = s(coset(ax
+1,ay
,az
),1,l
) -&
123 rab(1)*s(coset(ax
,ay
,az
),1,l
)
124 s(coset(ax
,ay
,az
),4,l
) = s(coset(ax
,ay
,az
+1),1,l
) -&
125 rab(3)*s(coset(ax
,ay
,az
),1,l
)
132 s(coset(ax
,ay
,az
),2,l
) = rbp(1)*s(coset(ax
,ay
,az
),1,l
)
134 s(coset(ax
,ay
,az
),2,l
) = rbp(1)*s(coset(ax
,ay
,az
),1,l
) +&
135 fx
*s(coset(ax
-1,ay
,az
),1,l
)
137 IF (lx1
> 0) s(coset(ax
,ay
,az
),2,l
) = s(coset(ax
,ay
,az
),2,l
) +&
138 f2x
*s(coset(ax
,ay
,az
),1,lx1
)
140 s(coset(ax
,ay
,az
),3,l
) = rbp(2)*s(coset(ax
,ay
,az
),1,l
)
142 s(coset(ax
,ay
,az
),3,l
) = rbp(2)*s(coset(ax
,ay
,az
),1,l
) +&
143 fy
*s(coset(ax
,ay
-1,az
),1,l
)
145 IF (ly1
> 0) s(coset(ax
,ay
,az
),3,l
) = s(coset(ax
,ay
,az
),3,l
) +&
146 f2y
*s(coset(ax
,ay
,az
),1,ly1
)
148 s(coset(ax
,ay
,az
),4,l
) = rbp(3)*s(coset(ax
,ay
,az
),1,l
)
150 s(coset(ax
,ay
,az
),4,l
) = rbp(3)*s(coset(ax
,ay
,az
),1,l
) +&
151 fz
*s(coset(ax
,ay
,az
-1),1,l
)
153 IF (lz1
> 0) s(coset(ax
,ay
,az
),4,l
) = s(coset(ax
,ay
,az
),4,l
) +&
154 f2z
*s(coset(ax
,ay
,az
),1,lz1
)
158 DO la
=la_start
,la_max
-1
161 s(coset(ax
,ay
,az
),coset(0,0,lb
),l
) =&
162 rab(3)*s(coset(ax
,ay
,az
),coset(0,0,lb
-1),l
)
165 s(coset(ax
,ay
,az
),coset(bx
,by
,bz
),l
) =&
166 rab(1)*s(coset(ax
,ay
,az
),coset(bx
-1,by
,bz
),l
)
175 s(coset(ax
,ay
,az
),coset(0,0,lb
),l
) =&
176 rbp(3)*s(coset(ax
,ay
,az
),coset(0,0,lb
-1),l
) +&
177 f3
*s(coset(ax
,ay
,az
),coset(0,0,lb
-2),l
)
179 IF (lz1
> 0) s(coset(ax
,ay
,az
),coset(0,0,lb
),l
) =&
180 f2z
*s(coset(ax
,ay
,az
),coset(0,0,lb
-1),lz1
)
182 IF (ly1
> 0) s(coset(ax
,ay
,az
),coset(0,1,bz
),l
) =&
183 f2y
*s(coset(ax
,ay
,az
),coset(0,0,bz
),ly1
)
185 s(coset(ax
,ay
,az
),coset(0,by
,bz
),l
) =&
186 f3
*s(coset(ax
,ay
,az
),coset(0,by
-2,bz
),l
)
187 IF (ly1
> 0) s(coset(ax
,ay
,az
),coset(0,by
,bz
),l
) =&
188 f2y
*s(coset(ax
,ay
,az
),coset(0,by
-1,bz
),ly1
)
190 s(coset(ax
,ay
,az
),coset(0,1,bz
),l
) =&
191 fy
*s(coset(ax
,ay
-1,az
),coset(0,0,bz
),l
)
195 IF (lx1
> 0) s(coset(ax
,ay
,az
),coset(1,by
,bz
),l
) =&
196 f2x
*s(coset(ax
,ay
,az
),coset(0,by
,bz
),lx1
)
200 s(coset(ax
,ay
,az
),coset(bx
,by
,bz
),l
) =&
201 f3
*s(coset(ax
,ay
,az
),coset(bx
-2,by
,bz
),l
)
202 IF (lx1
> 0) s(coset(ax
,ay
,az
),coset(bx
,by
,bz
),l
) =&
203 f2x
*s(coset(ax
,ay
,az
),coset(bx
-1,by
,bz
),lx1
)
207 IF (lx1
> 0) s(coset(ax
,ay
,az
),coset(1,by
,bz
),l
) =&
208 f2x
*s(coset(ax
,ay
,az
),coset(0,by
,bz
),lx1
)
212 s(coset(ax
,ay
,az
),coset(bx
,by
,bz
),l
) =&
213 f3
*s(coset(ax
,ay
,az
),coset(bx
-2,by
,bz
),l
)
214 IF (lx1
> 0) s(coset(ax
,ay
,az
),coset(bx
,by
,bz
),l
) =&
215 f2x
*s(coset(ax
,ay
,az
),coset(bx
-1,by
,bz
),lx1
)
225 IF (lz1
> 0) s(1,coset(0,0,lb
),l
) = s(1,coset(0,0,lb
),l
) +&
226 f2z
*s(1,coset(0,0,lb
-1),lz1
)
227 IF (ly1
> 0) s(1,coset(0,1,bz
),l
) = s(1,coset(0,1,bz
),l
) +&
228 f2y
*s(1,coset(0,0,bz
),ly1
)
230 s(1,coset(0,by
,bz
),l
) = rbp(2)*s(1,coset(0,by
-1,bz
),l
) +&
231 f2
*REAL(by
-1,dp
)*s(1,coset(0,by
-2,bz
),l
)
232 IF (lx1
> 0) s(1,coset(1,by
,bz
),l
) = s(1,coset(1,by
,bz
),l
) +&
233 f2x
*s(1,coset(0,by
,bz
),lx1
)
237 IF (lx1
> 0) s(1,coset(bx
,by
,bz
),l
) = s(1,coset(bx
,by
,bz
),l
) +&
238 f2x
*s(1,coset(bx
-1,by
,bz
),lx1
)
245 DO k
=2,ncoset(lc_max
)
246 DO j
=1,ncoset(lb_max
)
251 END SUBROUTINE moment
252 SUBROUTINE diff_momop(la_max
,npgfa
,zeta
,rpgfa
,la_min
,&
253 order
,rac
,rbc
,difmab
,mab_ext
)
254 REAL(KIND
=dp
), DIMENSION(:, :, :), &
255 OPTIONAL
, POINTER :: mab_ext
256 REAL(KIND
=dp
), ALLOCATABLE
, &
257 DIMENSION(:, :, :) :: difmab_tmp
258 DO imom
= 1,ncoset(order
)-1
259 CALL adbdr(la_max
,npgfa
,rpgfa
,la_min
,&
260 difmab_tmp(:,:,2), difmab_tmp(:,:,3))
262 END SUBROUTINE diff_momop
263 END MODULE ai_moments