1! { dg-do compile } 2! { dg-options "-O -fbounds-check -w" } 3MODULE kinds 4 INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND ( 14, 200 ) 5 INTEGER, DIMENSION(:), ALLOCATABLE :: nco,ncoset,nso,nsoset 6 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: co,coset 7END MODULE kinds 8MODULE ai_moments 9 USE kinds 10CONTAINS 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 16 DO ipgf=1,npgfa 17 DO jpgf=1,npgfb 18 IF (la_max > 0) THEN 19 DO la=2,la_max 20 DO ax=2,la 21 DO ay=0,la-ax 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) 26 END DO 27 END DO 28 END DO 29 IF (lb_max > 0) THEN 30 DO lb=2,lb_max 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)) 33 DO bx=2,lb 34 DO by=0,lb-bx 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)) 37 END DO 38 END DO 39 END DO 40 END IF 41 END IF 42 DO j=ncoset(lb_min-1)+1,ncoset(lb_max) 43 END DO 44 END DO 45 END DO 46 END SUBROUTINE cossin 47 SUBROUTINE moment(la_max,npgfa,zeta,rpgfa,la_min,& 48 lb_max,npgfb,zetb,rpgfb,& 49 lc_max,rac,rbc,mab) 50 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa 51 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb 52 REAL(KIND=dp), DIMENSION(:, :, :), & 53 INTENT(INOUT) :: mab 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 57 DO ipgf=1,npgfa 58 DO jpgf=1,npgfb 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) 63 mab(i,j,k) = 0.0_dp 64 END DO 65 END DO 66 END DO 67 END IF 68 rpc = zetp*(zeta(ipgf)*rac+zetb(jpgf)*rbc) 69 DO l=2, ncoset(lc_max) 70 lx = indco(1,l) 71 l2 = 0 72 IF ( lz > 0 ) THEN 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) 77 END IF 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) 80 END DO 81 DO l = 1, ncoset(lc_max) 82 IF ( lx > 0 ) THEN 83 lx1 = coset(lx-1,ly,lz) 84 END IF 85 IF ( ly > 0 ) THEN 86 ly1 = coset(lx,ly-1,lz) 87 END IF 88 IF (la_max > 0) THEN 89 DO la=2,la_max 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) 94 DO ay=2,la 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) 99 END DO 100 DO ay=0,la-1 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) 103 END DO 104 DO ax=2,la 105 DO ay=0,la-ax 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) 110 END DO 111 END DO 112 END DO 113 IF (lb_max > 0) THEN 114 DO j=2,ncoset(lb_max) 115 DO i=1,ncoset(la_max) 116 s(i,j,l) = 0.0_dp 117 END DO 118 END DO 119 DO la=la_start,la_max-1 120 DO ax=0,la 121 DO ay=0,la-ax 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) 126 END DO 127 END DO 128 END DO 129 DO ax=0,la_max 130 DO ay=0,la_max-ax 131 IF (ax == 0) THEN 132 s(coset(ax,ay,az),2,l) = rbp(1)*s(coset(ax,ay,az),1,l) 133 ELSE 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) 136 END IF 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) 139 IF (ay == 0) THEN 140 s(coset(ax,ay,az),3,l) = rbp(2)*s(coset(ax,ay,az),1,l) 141 ELSE 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) 144 END IF 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) 147 IF (az == 0) THEN 148 s(coset(ax,ay,az),4,l) = rbp(3)*s(coset(ax,ay,az),1,l) 149 ELSE 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) 152 END IF 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) 155 END DO 156 END DO 157 DO lb=2,lb_max 158 DO la=la_start,la_max-1 159 DO ax=0,la 160 DO ay=0,la-ax 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) 163 DO bx=1,lb 164 DO by=0,lb-bx 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) 167 END DO 168 END DO 169 END DO 170 END DO 171 END DO 172 DO ax=0,la_max 173 DO ay=0,la_max-ax 174 IF (az == 0) THEN 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) 178 END IF 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) 181 IF (ay == 0) THEN 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) 184 DO by=2,lb 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) 189 END DO 190 s(coset(ax,ay,az),coset(0,1,bz),l) =& 191 fy*s(coset(ax,ay-1,az),coset(0,0,bz),l) 192 END IF 193 IF (ax == 0) THEN 194 DO by=0,lb-1 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) 197 END DO 198 DO bx=2,lb 199 DO by=0,lb-bx 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) 204 END DO 205 END DO 206 DO by=0,lb-1 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) 209 END DO 210 DO bx=2,lb 211 DO by=0,lb-bx 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) 216 END DO 217 END DO 218 END IF 219 END DO 220 END DO 221 END DO 222 END IF 223 IF (lb_max > 0) THEN 224 DO lb=2,lb_max 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) 229 DO by=2,lb 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) 234 END DO 235 DO bx=2,lb 236 DO by=0,lb-bx 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) 239 END DO 240 END DO 241 END DO 242 END IF 243 END IF 244 END DO 245 DO k=2,ncoset(lc_max) 246 DO j=1,ncoset(lb_max) 247 END DO 248 END DO 249 END DO 250 END DO 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)) 261 END DO 262 END SUBROUTINE diff_momop 263END MODULE ai_moments 264