78 SUBROUTINE sptgpsv(IROMB,MAXWV,KMAX,NPS,
79 & KWSKIP,KGSKIP,NISKIP,NJSKIP,
80 & TRUE,XMESH,ORIENT,WAVED,WAVEZ,UN,VN,US,VS)
82 REAL WAVED(*),WAVEZ(*),UN(*),VN(*),US(*),VS(*)
83 REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
84 REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
85 REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
86 REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
88 REAL SLON(MAXWV,8),CLON(MAXWV,8),SROT(0:3),CROT(0:3)
89 REAL W((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2+1,2*KMAX)
90 REAL WTOP(2*(MAXWV+1),2*KMAX)
91 REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),PLNTOP(MAXWV+1)
92 REAL F(2*MAXWV+3,2,2*KMAX)
93 DATA srot/0.,1.,0.,-1./,crot/1.,0.,-1.,0./
94 parameter(rerth=6.3712e6)
95 parameter(pi=3.14159265358979,dpr=180./pi)
98 CALL spwget(iromb,maxwv,eps,epstop,enn1,elonn1,eon,eontop)
99 mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
108 IF(kg.EQ.0) kg=nps*nps
113 gq=((1.+sin(true/dpr))*rerth/xmesh)**2
120 CALL spdz2uv(iromb,maxwv,enn1,elonn1,eon,eontop,
121 & waved(kws+1),wavez(kws+1),
122 & w(1,k),w(1,kmax+k),wtop(1,k),wtop(1,kmax+k))
128 ij1=(i1-1)*ni+(j1-1)*nj+1
131 CALL splegend(iromb,maxwv,slat1,clat1,eps,epstop,
133 CALL spsynth(iromb,maxwv,2*maxwv,idim,mdim,2*mxtop,2*kmax,
134 & clat1,pln,plntop,mp,w,wtop,f)
142 un(ijk1)=2*( coso*f(3,1,ku)+sino*f(3,1,kv))
143 vn(ijk1)=2*(-sino*f(3,1,ku)+coso*f(3,1,kv))
144 us(ijk1)=2*( coso*f(3,2,ku)-sino*f(3,2,kv))
145 vs(ijk1)=2*( sino*f(3,2,ku)+coso*f(3,2,kv))
164 ij1=(i1-1)*ni+(j1-1)*nj+1
165 ij3=(i3-1)*ni+(j3-1)*nj+1
166 ij5=(i5-1)*ni+(j5-1)*nj+1
167 ij7=(i7-1)*ni+(j7-1)*nj+1
171 slat1=(gq-rq)/(gq+rq)
172 clat1=sqrt(1.-slat1**2)
173 CALL splegend(iromb,maxwv,slat1,clat1,eps,epstop,
175 CALL spsynth(iromb,maxwv,2*maxwv,idim,mdim,2*mxtop,2*kmax,
176 & clat1,pln,plntop,mp,w,wtop,f)
178 slon(l,1)=sin(l*radlon)
179 clon(l,1)=cos(l*radlon)
180 slon(l,3)=slon(l,1)*crot(mod(1*l,4))
181 & -clon(l,1)*srot(mod(1*l,4))
182 clon(l,3)=clon(l,1)*crot(mod(1*l,4))
183 & +slon(l,1)*srot(mod(1*l,4))
184 slon(l,5)=slon(l,1)*crot(mod(2*l,4))
185 & -clon(l,1)*srot(mod(2*l,4))
186 clon(l,5)=clon(l,1)*crot(mod(2*l,4))
187 & +slon(l,1)*srot(mod(2*l,4))
188 slon(l,7)=slon(l,1)*crot(mod(3*l,4))
189 & -clon(l,1)*srot(mod(3*l,4))
190 clon(l,7)=clon(l,1)*crot(mod(3*l,4))
191 & +slon(l,1)*srot(mod(3*l,4))
224 un(ij1)=un(ij1)+2*(f(lr,1,ku)*clon(l,1)
225 & -f(li,1,ku)*slon(l,1))
226 vn(ij1)=vn(ij1)+2*(f(lr,1,kv)*clon(l,1)
227 & -f(li,1,kv)*slon(l,1))
228 un(ij3)=un(ij3)+2*(f(lr,1,kv)*clon(l,3)
229 & -f(li,1,kv)*slon(l,3))
230 vn(ij3)=vn(ij3)-2*(f(lr,1,ku)*clon(l,3)
231 & -f(li,1,ku)*slon(l,3))
232 un(ij5)=un(ij5)-2*(f(lr,1,ku)*clon(l,5)
233 & -f(li,1,ku)*slon(l,5))
234 vn(ij5)=vn(ij5)-2*(f(lr,1,kv)*clon(l,5)
235 & -f(li,1,kv)*slon(l,5))
236 un(ij7)=un(ij7)-2*(f(lr,1,kv)*clon(l,7)
237 & -f(li,1,kv)*slon(l,7))
238 vn(ij7)=vn(ij7)+2*(f(lr,1,ku)*clon(l,7)
239 & -f(li,1,ku)*slon(l,7))
240 us(ij1)=us(ij1)-2*(f(lr,2,ku)*clon(l,5)
241 & -f(li,2,ku)*slon(l,5))
242 vs(ij1)=vs(ij1)-2*(f(lr,2,kv)*clon(l,5)
243 & -f(li,2,kv)*slon(l,5))
244 us(ij3)=us(ij3)-2*(f(lr,2,kv)*clon(l,3)
245 & -f(li,2,kv)*slon(l,3))
246 vs(ij3)=vs(ij3)+2*(f(lr,2,ku)*clon(l,3)
247 & -f(li,2,ku)*slon(l,3))
248 us(ij5)=us(ij5)+2*(f(lr,2,ku)*clon(l,1)
249 & -f(li,2,ku)*slon(l,1))
250 vs(ij5)=vs(ij5)+2*(f(lr,2,kv)*clon(l,1)
251 & -f(li,2,kv)*slon(l,1))
252 us(ij7)=us(ij7)+2*(f(lr,2,kv)*clon(l,7)
253 & -f(li,2,kv)*slon(l,7))
254 vs(ij7)=vs(ij7)-2*(f(lr,2,ku)*clon(l,7)
255 & -f(li,2,ku)*slon(l,7))
269 un(ijk1)=un(ijk1)+2*(f(lr,1,ku)*clon(l,1)
270 & -f(li,1,ku)*slon(l,1))
271 vn(ijk1)=vn(ijk1)+2*(f(lr,1,kv)*clon(l,1)
272 & -f(li,1,kv)*slon(l,1))
273 un(ijk3)=un(ijk3)+2*(f(lr,1,kv)*clon(l,3)
274 & -f(li,1,kv)*slon(l,3))
275 vn(ijk3)=vn(ijk3)-2*(f(lr,1,ku)*clon(l,3)
276 & -f(li,1,ku)*slon(l,3))
277 un(ijk5)=un(ijk5)-2*(f(lr,1,ku)*clon(l,5)
278 & -f(li,1,ku)*slon(l,5))
279 vn(ijk5)=vn(ijk5)-2*(f(lr,1,kv)*clon(l,5)
280 & -f(li,1,kv)*slon(l,5))
281 un(ijk7)=un(ijk7)-2*(f(lr,1,kv)*clon(l,7)
282 & -f(li,1,kv)*slon(l,7))
283 vn(ijk7)=vn(ijk7)+2*(f(lr,1,ku)*clon(l,7)
284 & -f(li,1,ku)*slon(l,7))
285 us(ijk1)=us(ijk1)-2*(f(lr,2,ku)*clon(l,5)
286 & -f(li,2,ku)*slon(l,5))
287 vs(ijk1)=vs(ijk1)-2*(f(lr,2,kv)*clon(l,5)
288 & -f(li,2,kv)*slon(l,5))
289 us(ijk3)=us(ijk3)-2*(f(lr,2,kv)*clon(l,3)
290 & -f(li,2,kv)*slon(l,3))
291 vs(ijk3)=vs(ijk3)+2*(f(lr,2,ku)*clon(l,3)
292 & -f(li,2,ku)*slon(l,3))
293 us(ijk5)=us(ijk5)+2*(f(lr,2,ku)*clon(l,1)
294 & -f(li,2,ku)*slon(l,1))
295 vs(ijk5)=vs(ijk5)+2*(f(lr,2,kv)*clon(l,1)
296 & -f(li,2,kv)*slon(l,1))
297 us(ijk7)=us(ijk7)+2*(f(lr,2,kv)*clon(l,7)
298 & -f(li,2,kv)*slon(l,7))
299 vs(ijk7)=vs(ijk7)-2*(f(lr,2,ku)*clon(l,7)
300 & -f(li,2,ku)*slon(l,7))
315 radlon=(orient-45)/dpr
322 ij1=(i1-1)*ni+(j1-1)*nj+1
323 ij3=(i3-1)*ni+(j3-1)*nj+1
324 ij5=(i5-1)*ni+(j5-1)*nj+1
325 ij7=(i7-1)*ni+(j7-1)*nj+1
329 slat1=(gq-rq)/(gq+rq)
330 clat1=sqrt(1.-slat1**2)
331 CALL splegend(iromb,maxwv,slat1,clat1,eps,epstop,
333 CALL spsynth(iromb,maxwv,2*maxwv,idim,mdim,2*mxtop,2*kmax,
334 & clat1,pln,plntop,mp,w,wtop,f)
336 slon(l,1)=sin(l*radlon)
337 clon(l,1)=cos(l*radlon)
338 slon(l,3)=slon(l,1)*crot(mod(1*l,4))
339 & -clon(l,1)*srot(mod(1*l,4))
340 clon(l,3)=clon(l,1)*crot(mod(1*l,4))
341 & +slon(l,1)*srot(mod(1*l,4))
342 slon(l,5)=slon(l,1)*crot(mod(2*l,4))
343 & -clon(l,1)*srot(mod(2*l,4))
344 clon(l,5)=clon(l,1)*crot(mod(2*l,4))
345 & +slon(l,1)*srot(mod(2*l,4))
346 slon(l,7)=slon(l,1)*crot(mod(3*l,4))
347 & -clon(l,1)*srot(mod(3*l,4))
348 clon(l,7)=clon(l,1)*crot(mod(3*l,4))
349 & +slon(l,1)*srot(mod(3*l,4))
359 un(ijk1)=srh*( f(1,1,ku)+f(1,1,kv))
360 vn(ijk1)=srh*(-f(1,1,ku)+f(1,1,kv))
361 un(ijk3)=srh*(-f(1,1,ku)+f(1,1,kv))
362 vn(ijk3)=srh*(-f(1,1,ku)-f(1,1,kv))
363 un(ijk5)=srh*(-f(1,1,ku)-f(1,1,kv))
364 vn(ijk5)=srh*( f(1,1,ku)-f(1,1,kv))
365 un(ijk7)=srh*( f(1,1,ku)-f(1,1,kv))
366 vn(ijk7)=srh*( f(1,1,ku)+f(1,1,kv))
367 us(ijk1)=srh*(-f(1,2,ku)-f(1,2,kv))
368 vs(ijk1)=srh*( f(1,2,ku)-f(1,2,kv))
369 us(ijk3)=srh*( f(1,2,ku)-f(1,2,kv))
370 vs(ijk3)=srh*( f(1,2,ku)+f(1,2,kv))
371 us(ijk5)=srh*( f(1,2,ku)+f(1,2,kv))
372 vs(ijk5)=srh*(-f(1,2,ku)+f(1,2,kv))
373 us(ijk7)=srh*(-f(1,2,ku)+f(1,2,kv))
374 vs(ijk7)=srh*(-f(1,2,ku)-f(1,2,kv))
382 un(ij1)=un(ij1)+2*srh*(( f(lr,1,ku)+f(lr,1,kv))
384 & -( f(li,1,ku)+f(li,1,kv))
386 vn(ij1)=vn(ij1)+2*srh*((-f(lr,1,ku)+f(lr,1,kv))
388 & -(-f(li,1,ku)+f(li,1,kv))
390 un(ij3)=un(ij3)+2*srh*((-f(lr,1,ku)+f(lr,1,kv))
392 & -(-f(li,1,ku)+f(li,1,kv))
394 vn(ij3)=vn(ij3)+2*srh*((-f(lr,1,ku)-f(lr,1,kv))
396 & -(-f(li,1,ku)-f(li,1,kv))
398 un(ij5)=un(ij5)+2*srh*((-f(lr,1,ku)-f(lr,1,kv))
400 & -(-f(li,1,ku)-f(li,1,kv))
402 vn(ij5)=vn(ij5)+2*srh*(( f(lr,1,ku)-f(lr,1,kv))
404 & -( f(li,1,ku)-f(li,1,kv))
406 un(ij7)=un(ij7)+2*srh*(( f(lr,1,ku)-f(lr,1,kv))
408 & -( f(li,1,ku)-f(li,1,kv))
410 vn(ij7)=vn(ij7)+2*srh*(( f(lr,1,ku)+f(lr,1,kv))
412 & -( f(li,1,ku)+f(li,1,kv))
414 us(ij1)=us(ij1)+2*srh*((-f(lr,2,ku)-f(lr,2,kv))
416 & -(-f(li,2,ku)-f(li,2,kv))
418 vs(ij1)=vs(ij1)+2*srh*(( f(lr,2,ku)-f(lr,2,kv))
420 & -( f(li,2,ku)-f(li,2,kv))
422 us(ij3)=us(ij3)+2*srh*(( f(lr,2,ku)-f(lr,2,kv))
424 & -( f(li,2,ku)-f(li,2,kv))
426 vs(ij3)=vs(ij3)+2*srh*(( f(lr,2,ku)+f(lr,2,kv))
428 & -( f(li,2,ku)+f(li,2,kv))
430 us(ij5)=us(ij5)+2*srh*(( f(lr,2,ku)+f(lr,2,kv))
432 & -( f(li,2,ku)+f(li,2,kv))
434 vs(ij5)=vs(ij5)+2*srh*((-f(lr,2,ku)+f(lr,2,kv))
436 & -(-f(li,2,ku)+f(li,2,kv))
438 us(ij7)=us(ij7)+2*srh*((-f(lr,2,ku)+f(lr,2,kv))
440 & -(-f(li,2,ku)+f(li,2,kv))
442 vs(ij7)=vs(ij7)+2*srh*((-f(lr,2,ku)-f(lr,2,kv))
444 & -(-f(li,2,ku)-f(li,2,kv))
459 un(ijk1)=un(ijk1)+2*srh*(( f(lr,1,ku)+f(lr,1,kv))
461 & -( f(li,1,ku)+f(li,1,kv))
463 vn(ijk1)=vn(ijk1)+2*srh*((-f(lr,1,ku)+f(lr,1,kv))
465 & -(-f(li,1,ku)+f(li,1,kv))
467 un(ijk3)=un(ijk3)+2*srh*((-f(lr,1,ku)+f(lr,1,kv))
469 & -(-f(li,1,ku)+f(li,1,kv))
471 vn(ijk3)=vn(ijk3)+2*srh*((-f(lr,1,ku)-f(lr,1,kv))
473 & -(-f(li,1,ku)-f(li,1,kv))
475 un(ijk5)=un(ijk5)+2*srh*((-f(lr,1,ku)-f(lr,1,kv))
477 & -(-f(li,1,ku)-f(li,1,kv))
479 vn(ijk5)=vn(ijk5)+2*srh*(( f(lr,1,ku)-f(lr,1,kv))
481 & -( f(li,1,ku)-f(li,1,kv))
483 un(ijk7)=un(ijk7)+2*srh*(( f(lr,1,ku)-f(lr,1,kv))
485 & -( f(li,1,ku)-f(li,1,kv))
487 vn(ijk7)=vn(ijk7)+2*srh*(( f(lr,1,ku)+f(lr,1,kv))
489 & -( f(li,1,ku)+f(li,1,kv))
491 us(ijk1)=us(ijk1)+2*srh*((-f(lr,2,ku)-f(lr,2,kv))
493 & -(-f(li,2,ku)-f(li,2,kv))
495 vs(ijk1)=vs(ijk1)+2*srh*(( f(lr,2,ku)-f(lr,2,kv))
497 & -( f(li,2,ku)-f(li,2,kv))
499 us(ijk3)=us(ijk3)+2*srh*(( f(lr,2,ku)-f(lr,2,kv))
501 & -( f(li,2,ku)-f(li,2,kv))
503 vs(ijk3)=vs(ijk3)+2*srh*(( f(lr,2,ku)+f(lr,2,kv))
505 & -( f(li,2,ku)+f(li,2,kv))
507 us(ijk5)=us(ijk5)+2*srh*(( f(lr,2,ku)+f(lr,2,kv))
509 & -( f(li,2,ku)+f(li,2,kv))
511 vs(ijk5)=vs(ijk5)+2*srh*((-f(lr,2,ku)+f(lr,2,kv))
513 & -(-f(li,2,ku)+f(li,2,kv))
515 us(ijk7)=us(ijk7)+2*srh*((-f(lr,2,ku)+f(lr,2,kv))
517 & -(-f(li,2,ku)+f(li,2,kv))
519 vs(ijk7)=vs(ijk7)+2*srh*((-f(lr,2,ku)-f(lr,2,kv))
521 & -(-f(li,2,ku)-f(li,2,kv))
552 ij1=(i1-1)*ni+(j1-1)*nj+1
553 ij2=(i2-1)*ni+(j2-1)*nj+1
554 ij3=(i3-1)*ni+(j3-1)*nj+1
555 ij4=(i4-1)*ni+(j4-1)*nj+1
556 ij5=(i5-1)*ni+(j5-1)*nj+1
557 ij6=(i6-1)*ni+(j6-1)*nj+1
558 ij7=(i7-1)*ni+(j7-1)*nj+1
559 ij8=(i8-1)*ni+(j8-1)*nj+1
564 slat1=(gq-rq)/(gq+rq)
565 clat1=sqrt(1.-slat1**2)
566 radlon1=orient/dpr+atan(-di1/dj1)
567 radlon2=(orient-45)/dpr*2-radlon1
568 CALL splegend(iromb,maxwv,slat1,clat1,eps,epstop,
570 CALL spsynth(iromb,maxwv,2*maxwv,idim,mdim,2*mxtop,2*kmax,
571 & clat1,pln,plntop,mp,w,wtop,f)
573 slon(l,1)=sin(l*radlon1)
574 clon(l,1)=cos(l*radlon1)
575 slon(l,2)=sin(l*radlon2)
576 clon(l,2)=cos(l*radlon2)
577 slon(l,3)=slon(l,1)*crot(mod(1*l,4))
578 & -clon(l,1)*srot(mod(1*l,4))
579 clon(l,3)=clon(l,1)*crot(mod(1*l,4))
580 & +slon(l,1)*srot(mod(1*l,4))
581 slon(l,4)=slon(l,2)*crot(mod(1*l,4))
582 & -clon(l,2)*srot(mod(1*l,4))
583 clon(l,4)=clon(l,2)*crot(mod(1*l,4))
584 & +slon(l,2)*srot(mod(1*l,4))
585 slon(l,5)=slon(l,1)*crot(mod(2*l,4))
586 & -clon(l,1)*srot(mod(2*l,4))
587 clon(l,5)=clon(l,1)*crot(mod(2*l,4))
588 & +slon(l,1)*srot(mod(2*l,4))
589 slon(l,6)=slon(l,2)*crot(mod(2*l,4))
590 & -clon(l,2)*srot(mod(2*l,4))
591 clon(l,6)=clon(l,2)*crot(mod(2*l,4))
592 & +slon(l,2)*srot(mod(2*l,4))
593 slon(l,7)=slon(l,1)*crot(mod(3*l,4))
594 & -clon(l,1)*srot(mod(3*l,4))
595 clon(l,7)=clon(l,1)*crot(mod(3*l,4))
596 & +slon(l,1)*srot(mod(3*l,4))
597 slon(l,8)=slon(l,2)*crot(mod(3*l,4))
598 & -clon(l,2)*srot(mod(3*l,4))
599 clon(l,8)=clon(l,2)*crot(mod(3*l,4))
600 & +slon(l,2)*srot(mod(3*l,4))
614 un(ijk1)=rr*(-dj1*f(1,1,ku)-di1*f(1,1,kv))
615 vn(ijk1)=rr*( di1*f(1,1,ku)-dj1*f(1,1,kv))
616 un(ijk2)=rr*(-di1*f(1,1,ku)-dj1*f(1,1,kv))
617 vn(ijk2)=rr*( dj1*f(1,1,ku)-di1*f(1,1,kv))
618 un(ijk3)=rr*( di1*f(1,1,ku)-dj1*f(1,1,kv))
619 vn(ijk3)=rr*( dj1*f(1,1,ku)+di1*f(1,1,kv))
620 un(ijk4)=rr*( dj1*f(1,1,ku)-di1*f(1,1,kv))
621 vn(ijk4)=rr*( di1*f(1,1,ku)+dj1*f(1,1,kv))
622 un(ijk5)=rr*( dj1*f(1,1,ku)+di1*f(1,1,kv))
623 vn(ijk5)=rr*(-di1*f(1,1,ku)+dj1*f(1,1,kv))
624 un(ijk6)=rr*( di1*f(1,1,ku)+dj1*f(1,1,kv))
625 vn(ijk6)=rr*(-dj1*f(1,1,ku)+di1*f(1,1,kv))
626 un(ijk7)=rr*(-di1*f(1,1,ku)+dj1*f(1,1,kv))
627 vn(ijk7)=rr*(-dj1*f(1,1,ku)-di1*f(1,1,kv))
628 un(ijk8)=rr*(-dj1*f(1,1,ku)+di1*f(1,1,kv))
629 vn(ijk8)=rr*(-di1*f(1,1,ku)-dj1*f(1,1,kv))
630 us(ijk1)=rr*( dj1*f(1,2,ku)+di1*f(1,2,kv))
631 vs(ijk1)=rr*(-di1*f(1,2,ku)+dj1*f(1,2,kv))
632 us(ijk2)=rr*( di1*f(1,2,ku)+dj1*f(1,2,kv))
633 vs(ijk2)=rr*(-dj1*f(1,2,ku)+di1*f(1,2,kv))
634 us(ijk3)=rr*(-di1*f(1,2,ku)+dj1*f(1,2,kv))
635 vs(ijk3)=rr*(-dj1*f(1,2,ku)-di1*f(1,2,kv))
636 us(ijk4)=rr*(-dj1*f(1,2,ku)+di1*f(1,2,kv))
637 vs(ijk4)=rr*(-di1*f(1,2,ku)-dj1*f(1,2,kv))
638 us(ijk5)=rr*(-dj1*f(1,2,ku)-di1*f(1,2,kv))
639 vs(ijk5)=rr*( di1*f(1,2,ku)-dj1*f(1,2,kv))
640 us(ijk6)=rr*(-di1*f(1,2,ku)-dj1*f(1,2,kv))
641 vs(ijk6)=rr*( dj1*f(1,2,ku)-di1*f(1,2,kv))
642 us(ijk7)=rr*( di1*f(1,2,ku)-dj1*f(1,2,kv))
643 vs(ijk7)=rr*( dj1*f(1,2,ku)+di1*f(1,2,kv))
644 us(ijk8)=rr*( dj1*f(1,2,ku)-di1*f(1,2,kv))
645 vs(ijk8)=rr*( di1*f(1,2,ku)+dj1*f(1,2,kv))
653 un(ij1)=un(ij1)+2*rr*((-dj1*f(lr,1,ku)-di1*f(lr,1,kv))
655 & -(-dj1*f(li,1,ku)-di1*f(li,1,kv))
657 vn(ij1)=vn(ij1)+2*rr*(( di1*f(lr,1,ku)-dj1*f(lr,1,kv))
659 & -( di1*f(li,1,ku)-dj1*f(li,1,kv))
661 un(ij2)=un(ij2)+2*rr*((-di1*f(lr,1,ku)-dj1*f(lr,1,kv))
663 & -(-di1*f(li,1,ku)-dj1*f(li,1,kv))
665 vn(ij2)=vn(ij2)+2*rr*(( dj1*f(lr,1,ku)-di1*f(lr,1,kv))
667 & -( dj1*f(li,1,ku)-di1*f(li,1,kv))
669 un(ij3)=un(ij3)+2*rr*(( di1*f(lr,1,ku)-dj1*f(lr,1,kv))
671 & -( di1*f(li,1,ku)-dj1*f(li,1,kv))
673 vn(ij3)=vn(ij3)+2*rr*(( dj1*f(lr,1,ku)+di1*f(lr,1,kv))
675 & -( dj1*f(li,1,ku)+di1*f(li,1,kv))
677 un(ij4)=un(ij4)+2*rr*(( dj1*f(lr,1,ku)-di1*f(lr,1,kv))
679 & -( dj1*f(li,1,ku)-di1*f(li,1,kv))
681 vn(ij4)=vn(ij4)+2*rr*(( di1*f(lr,1,ku)+dj1*f(lr,1,kv))
683 & -( di1*f(li,1,ku)+dj1*f(li,1,kv))
685 un(ij5)=un(ij5)+2*rr*(( dj1*f(lr,1,ku)+di1*f(lr,1,kv))
687 & -( dj1*f(li,1,ku)+di1*f(li,1,kv))
689 vn(ij5)=vn(ij5)+2*rr*((-di1*f(lr,1,ku)+dj1*f(lr,1,kv))
691 & -(-di1*f(li,1,ku)+dj1*f(li,1,kv))
693 un(ij6)=un(ij6)+2*rr*(( di1*f(lr,1,ku)+dj1*f(lr,1,kv))
695 & -( di1*f(li,1,ku)+dj1*f(li,1,kv))
697 vn(ij6)=vn(ij6)+2*rr*((-dj1*f(lr,1,ku)+di1*f(lr,1,kv))
699 & -(-dj1*f(li,1,ku)+di1*f(li,1,kv))
701 un(ij7)=un(ij7)+2*rr*((-di1*f(lr,1,ku)+dj1*f(lr,1,kv))
703 & -(-di1*f(li,1,ku)+dj1*f(li,1,kv))
705 vn(ij7)=vn(ij7)+2*rr*((-dj1*f(lr,1,ku)-di1*f(lr,1,kv))
707 & -(-dj1*f(li,1,ku)-di1*f(li,1,kv))
709 un(ij8)=un(ij8)+2*rr*((-dj1*f(lr,1,ku)+di1*f(lr,1,kv))
711 & -(-dj1*f(li,1,ku)+di1*f(li,1,kv))
713 vn(ij8)=vn(ij8)+2*rr*((-di1*f(lr,1,ku)-dj1*f(lr,1,kv))
715 & -(-di1*f(li,1,ku)-dj1*f(li,1,kv))
717 us(ij1)=us(ij1)+2*rr*(( dj1*f(lr,2,ku)+di1*f(lr,2,kv))
719 & -( dj1*f(li,2,ku)+di1*f(li,2,kv))
721 vs(ij1)=vs(ij1)+2*rr*((-di1*f(lr,2,ku)+dj1*f(lr,2,kv))
723 & -(-di1*f(li,2,ku)+dj1*f(li,2,kv))
725 us(ij2)=us(ij2)+2*rr*(( di1*f(lr,2,ku)+dj1*f(lr,2,kv))
727 & -( di1*f(li,2,ku)+dj1*f(li,2,kv))
729 vs(ij2)=vs(ij2)+2*rr*((-dj1*f(lr,2,ku)+di1*f(lr,2,kv))
731 & -(-dj1*f(li,2,ku)+di1*f(li,2,kv))
733 us(ij3)=us(ij3)+2*rr*((-di1*f(lr,2,ku)+dj1*f(lr,2,kv))
735 & -(-di1*f(li,2,ku)+dj1*f(li,2,kv))
737 vs(ij3)=vs(ij3)+2*rr*((-dj1*f(lr,2,ku)-di1*f(lr,2,kv))
739 & -(-dj1*f(li,2,ku)-di1*f(li,2,kv))
741 us(ij4)=us(ij4)+2*rr*((-dj1*f(lr,2,ku)+di1*f(lr,2,kv))
743 & -(-dj1*f(li,2,ku)+di1*f(li,2,kv))
745 vs(ij4)=vs(ij4)+2*rr*((-di1*f(lr,2,ku)-dj1*f(lr,2,kv))
747 & -(-di1*f(li,2,ku)-dj1*f(li,2,kv))
749 us(ij5)=us(ij5)+2*rr*((-dj1*f(lr,2,ku)-di1*f(lr,2,kv))
751 & -(-dj1*f(li,2,ku)-di1*f(li,2,kv))
753 vs(ij5)=vs(ij5)+2*rr*(( di1*f(lr,2,ku)-dj1*f(lr,2,kv))
755 & -( di1*f(li,2,ku)-dj1*f(li,2,kv))
757 us(ij6)=us(ij6)+2*rr*((-di1*f(lr,2,ku)-dj1*f(lr,2,kv))
759 & -(-di1*f(li,2,ku)-dj1*f(li,2,kv))
761 vs(ij6)=vs(ij6)+2*rr*(( dj1*f(lr,2,ku)-di1*f(lr,2,kv))
763 & -( dj1*f(li,2,ku)-di1*f(li,2,kv))
765 us(ij7)=us(ij7)+2*rr*(( di1*f(lr,2,ku)-dj1*f(lr,2,kv))
767 & -( di1*f(li,2,ku)-dj1*f(li,2,kv))
769 vs(ij7)=vs(ij7)+2*rr*(( dj1*f(lr,2,ku)+di1*f(lr,2,kv))
771 & -( dj1*f(li,2,ku)+di1*f(li,2,kv))
773 us(ij8)=us(ij8)+2*rr*(( dj1*f(lr,2,ku)-di1*f(lr,2,kv))
775 & -( dj1*f(li,2,ku)-di1*f(li,2,kv))
777 vs(ij8)=vs(ij8)+2*rr*(( di1*f(lr,2,ku)+dj1*f(lr,2,kv))
779 & -( di1*f(li,2,ku)+dj1*f(li,2,kv))
798 un(ijk1)=un(ijk1)+2*rr*((-dj1*f(lr,1,ku)-di1*f(lr,1,kv))
800 & -(-dj1*f(li,1,ku)-di1*f(li,1,kv))
802 vn(ijk1)=vn(ijk1)+2*rr*(( di1*f(lr,1,ku)-dj1*f(lr,1,kv))
804 & -( di1*f(li,1,ku)-dj1*f(li,1,kv))
806 un(ijk2)=un(ijk2)+2*rr*((-di1*f(lr,1,ku)-dj1*f(lr,1,kv))
808 & -(-di1*f(li,1,ku)-dj1*f(li,1,kv))
810 vn(ijk2)=vn(ijk2)+2*rr*(( dj1*f(lr,1,ku)-di1*f(lr,1,kv))
812 & -( dj1*f(li,1,ku)-di1*f(li,1,kv))
814 un(ijk3)=un(ijk3)+2*rr*(( di1*f(lr,1,ku)-dj1*f(lr,1,kv))
816 & -( di1*f(li,1,ku)-dj1*f(li,1,kv))
818 vn(ijk3)=vn(ijk3)+2*rr*(( dj1*f(lr,1,ku)+di1*f(lr,1,kv))
820 & -( dj1*f(li,1,ku)+di1*f(li,1,kv))
822 un(ijk4)=un(ijk4)+2*rr*(( dj1*f(lr,1,ku)-di1*f(lr,1,kv))
824 & -( dj1*f(li,1,ku)-di1*f(li,1,kv))
826 vn(ijk4)=vn(ijk4)+2*rr*(( di1*f(lr,1,ku)+dj1*f(lr,1,kv))
828 & -( di1*f(li,1,ku)+dj1*f(li,1,kv))
830 un(ijk5)=un(ijk5)+2*rr*(( dj1*f(lr,1,ku)+di1*f(lr,1,kv))
832 & -( dj1*f(li,1,ku)+di1*f(li,1,kv))
834 vn(ijk5)=vn(ijk5)+2*rr*((-di1*f(lr,1,ku)+dj1*f(lr,1,kv))
836 & -(-di1*f(li,1,ku)+dj1*f(li,1,kv))
838 un(ijk6)=un(ijk6)+2*rr*(( di1*f(lr,1,ku)+dj1*f(lr,1,kv))
840 & -( di1*f(li,1,ku)+dj1*f(li,1,kv))
842 vn(ijk6)=vn(ijk6)+2*rr*((-dj1*f(lr,1,ku)+di1*f(lr,1,kv))
844 & -(-dj1*f(li,1,ku)+di1*f(li,1,kv))
846 un(ijk7)=un(ijk7)+2*rr*((-di1*f(lr,1,ku)+dj1*f(lr,1,kv))
848 & -(-di1*f(li,1,ku)+dj1*f(li,1,kv))
850 vn(ijk7)=vn(ijk7)+2*rr*((-dj1*f(lr,1,ku)-di1*f(lr,1,kv))
852 & -(-dj1*f(li,1,ku)-di1*f(li,1,kv))
854 un(ijk8)=un(ijk8)+2*rr*((-dj1*f(lr,1,ku)+di1*f(lr,1,kv))
856 & -(-dj1*f(li,1,ku)+di1*f(li,1,kv))
858 vn(ijk8)=vn(ijk8)+2*rr*((-di1*f(lr,1,ku)-dj1*f(lr,1,kv))
860 & -(-di1*f(li,1,ku)-dj1*f(li,1,kv))
862 us(ijk1)=us(ijk1)+2*rr*(( dj1*f(lr,2,ku)+di1*f(lr,2,kv))
864 & -( dj1*f(li,2,ku)+di1*f(li,2,kv))
866 vs(ijk1)=vs(ijk1)+2*rr*((-di1*f(lr,2,ku)+dj1*f(lr,2,kv))
868 & -(-di1*f(li,2,ku)+dj1*f(li,2,kv))
870 us(ijk2)=us(ijk2)+2*rr*(( di1*f(lr,2,ku)+dj1*f(lr,2,kv))
872 & -( di1*f(li,2,ku)+dj1*f(li,2,kv))
874 vs(ijk2)=vs(ijk2)+2*rr*((-dj1*f(lr,2,ku)+di1*f(lr,2,kv))
876 & -(-dj1*f(li,2,ku)+di1*f(li,2,kv))
878 us(ijk3)=us(ijk3)+2*rr*((-di1*f(lr,2,ku)+dj1*f(lr,2,kv))
880 & -(-di1*f(li,2,ku)+dj1*f(li,2,kv))
882 vs(ijk3)=vs(ijk3)+2*rr*((-dj1*f(lr,2,ku)-di1*f(lr,2,kv))
884 & -(-dj1*f(li,2,ku)-di1*f(li,2,kv))
886 us(ijk4)=us(ijk4)+2*rr*((-dj1*f(lr,2,ku)+di1*f(lr,2,kv))
888 & -(-dj1*f(li,2,ku)+di1*f(li,2,kv))
890 vs(ijk4)=vs(ijk4)+2*rr*((-di1*f(lr,2,ku)-dj1*f(lr,2,kv))
892 & -(-di1*f(li,2,ku)-dj1*f(li,2,kv))
894 us(ijk5)=us(ijk5)+2*rr*((-dj1*f(lr,2,ku)-di1*f(lr,2,kv))
896 & -(-dj1*f(li,2,ku)-di1*f(li,2,kv))
898 vs(ijk5)=vs(ijk5)+2*rr*(( di1*f(lr,2,ku)-dj1*f(lr,2,kv))
900 & -( di1*f(li,2,ku)-dj1*f(li,2,kv))
902 us(ijk6)=us(ijk6)+2*rr*((-di1*f(lr,2,ku)-dj1*f(lr,2,kv))
904 & -(-di1*f(li,2,ku)-dj1*f(li,2,kv))
906 vs(ijk6)=vs(ijk6)+2*rr*(( dj1*f(lr,2,ku)-di1*f(lr,2,kv))
908 & -( dj1*f(li,2,ku)-di1*f(li,2,kv))
910 us(ijk7)=us(ijk7)+2*rr*(( di1*f(lr,2,ku)-dj1*f(lr,2,kv))
912 & -( di1*f(li,2,ku)-dj1*f(li,2,kv))
914 vs(ijk7)=vs(ijk7)+2*rr*(( dj1*f(lr,2,ku)+di1*f(lr,2,kv))
916 & -( dj1*f(li,2,ku)+di1*f(li,2,kv))
918 us(ijk8)=us(ijk8)+2*rr*(( dj1*f(lr,2,ku)-di1*f(lr,2,kv))
920 & -( dj1*f(li,2,ku)-di1*f(li,2,kv))
922 vs(ijk8)=vs(ijk8)+2*rr*(( di1*f(lr,2,ku)+dj1*f(lr,2,kv))
924 & -( di1*f(li,2,ku)+dj1*f(li,2,kv))