1 SUBROUTINE furplm(p,grsq,x,y,v,xp,yp,gsq,np,nx,ny,ntx,nw,nh,is,js)
25 INTEGER :: js, is, ntx, ny, nx, np, nw, nh, nt, nx2, ny2
26 INTEGER :: jst, jdec, kprt, j, i, k, in, iq, id, i1, jy,
27 . ix, igp, l, iyy, ixx, ig, il
28 REAL(rprec) :: xp(*),yp(*)
29 REAL(rprec) :: grsq(nw,nh),gsq(*)
30 REAL(rprec) :: p(nw,*),x(*),y(*)
31 REAL(rprec) :: a, b, c, f, yq, v, dx, dy, dp, d, e, ff, ws,
41 IF(np.ne.-1)
go to 1999
45 1999
IF(nt.gt.0)
go to 2004
68 13
IF((v-p(i,j))*(p(i+1,j)-v).gt.0.)
go to 30
69 IF(abs(p(i,j)-v).ge.d)
go to 32
77 IF(jdec.lt.0)
go to 2005
81 2005
IF(j.lt.1)
go to 10
84 30
IF(i.gt.(nx-2))
go to 74
87 CALL fit(2,x(i-1),x(i),x(i+1),x(i+2),p(i-1,j),p(i,j),p(i+1,j),
88 1 p(i+2,j),xp(1),v,yq)
90 CALL cubic(x(i-1),x(i),x(i+1),x(i+2),ws,a,b,c,f)
91 gsq(1)=a*grsq(i-1,j)+b*grsq(i,j)+c*grsq(i+1,j)+f*grsq(i+2,j)
94 CALL fit(2,x(1),x(2),x(3),x(4),p(1,j),p(2,j),p(3,j),p(4,j),
97 CALL cubic(x(1),x(2),x(3),x(4),ws,a,b,c,f)
98 gsq(1)=a*grsq(1,j)+b*grsq(2,j)+c*grsq(3,j)+f*grsq(4,j)
101 CALL fit(2,x(nx-3),x(nx-2),x(nx-1),x(nx),p(nx-3,j),p(nx-2,j),
102 1 p(nx-1,j),p(nx,j),xp(1),v,yq)
104 CALL cubic(x(nx-3),x(nx-2),x(nx-1),x(nx),ws,a,b,c,f)
105 gsq(1)=a*grsq(nx-3,j)+b*grsq(nx-2,j)+c*grsq(nx-1,j)+f*grsq(nx,j)
108 IF(xp(1).le.x(i)) xp(1)=x(i)+ff
109 IF(xp(1).ge.x(i+1)) xp(1)=x(i+1)-ff
114 IF(kprt.eq.1)
WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
116 102
FORMAT(i5,3e13.6,10i5)
125 4
IF(il.eq.1)
go to 21
129 IF(i+1.gt.nx)
go to 10
133 IF(abs(v-p(i+1,j+1)).ge.d)
go to 45
140 45
IF(abs(v-p(i+1,j)).ge.d)
go to 33
149 33
IF((v-p(i,j))*(p(i+1,j)-v).le.0.)
go to 50
154 50
IF((v-p(i,j+1))*(p(i+1,j+1)-v).le.0)
go to 51
159 51
IF((v-p(i+1,j))*(p(i+1,j+1)-v).le.0.)
go to 52
168 IF(i-1.lt.0)
go to 10
172 21
IF(abs(v-p(i-1,j)).ge.d)
go to 46
178 46
IF(abs(v-p(i-1,j+1)).ge.d)
go to 34
188 34
IF((v-p(i,j+1))*(p(i-1,j+1)-v).le.0.)
go to 53
193 53
IF((v-p(i-1,j+1))*(p(i-1,j)-v).le.0)
go to 54
198 54
IF((v-p(i,j))*(p(i-1,j)-v).le.0.)
go to 55
210 5
IF(id.eq.15.and.il.eq.1) i=i-1
212 62
IF(in.eq.1)
go to 12
216 IF(j-1.lt.1)
go to 10
217 11
IF(abs(v-p(i+1,j-1)).ge.d)
go to 42
224 42
IF(abs(v-p(i,j-1)).ge.d)
go to 17
233 17
IF((v-p(i+1,j-1))*(p(i+1,j)-v).le.0.)
go to 56
234 IF(abs(v-p(i+1,j)).le.d)
go to 56
235 IF(id.eq.15.and.igp.eq.7)
go to 56
240 56
IF((v-p(i,j-1))*(p(i+1,j-1)-v).le.0.)
go to 57
245 57
IF((v-p(i,j))*(p(i,j-1)-v).le.0.)
go to 44
246 IF(abs(v-p(i,j)).le.d)
go to 44
251 44
IF(iq.ne.1)
go to 48
262 12
IF(j+1.gt.ny)
go to 10
267 IF(abs(v-p(i+1,j+1)).ge.d)
go to 40
274 40
IF(abs(v-p(i,j+1)).ge.d)
go to 61
283 61
IF((v-p(i,j+1))*(p(i+1,j+1)-v).le.0)
go to 59
288 59
IF((v-p(i+1,j))*(p(i+1,j+1)-v).le.0.)
go to 60
289 IF(abs(v-p(i+1,j)).le.d)
go to 60
294 60
IF((v-p(i,j+1))*(p(i,j)-v).le.0.)
go to 41
295 IF(abs(v-p(i,j)).le.d)
go to 41
300 41
IF(iq.ne.1)
go to 49
314 IF(k.eq.nt)
WRITE(66,106)
315 106
FORMAT(/,
"***nt points found in furplm***",/)
316 IF(k.eq.nt)
go to 2000
318 IF(ix.gt.(nx-2))
go to 3
320 CALL fit(2,x(ix-1),x(ix),x(ix+1),x(ix+2),p(ix-1,jy),p(ix,jy),
321 1 p(ix+1,jy),p(ix+2,jy),xp(k),v,yq)
323 CALL cubic(x(ix-1),x(ix),x(ix+1),x(ix+2),ws,a,b,c,f)
325 gsq(k)=a*grsq(ix-1,jy)+b*grsq(ix,jy)+c*grsq(ix+1,jy)
330 CALL fit(2,x(1),x(2),x(3),x(4),p(1,jy),p(2,jy),p(3,jy),p(4,jy),
333 CALL cubic(x(1),x(2),x(3),x(4),ws,a,b,c,f)
334 gsq(k)=a*grsq(1,jy)+b*grsq(2,jy)+c*grsq(3,jy)+f*grsq(4,jy)
337 CALL fit(2,x(nx-3),x(nx-2),x(nx-1),x(nx),p(nx-3,jy),p(nx-2,jy),
338 1 p(nx-1,jy),p(nx,jy),xp(k),v,yq)
340 CALL cubic(x(nx-3),x(nx-2),x(nx-1),x(nx),ws,a,b,c,f)
341 gsq(k)=a*grsq(nx-3,jy)+b*grsq(nx-2,jy)+c*grsq(nx-1,jy)
344 IF((id.eq.4).or.(id.eq.6)) i=i-1
345 IF(xp(k).le.x(ix)) xp(k)=x(ix)+ff
346 IF(xp(k).ge.x(ix+1)) xp(k)=x(ix+1)-ff
348 IF((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
349 1.and.k.ne.1)
go to 10
350 20
IF(jy.eq.1.and.(((id.eq.1).or.(id.eq.6)).or.(id.eq.8)))
go to 10
351 IF(jy.eq.ny.and.(((id.eq.2).or.(id.eq.4)).or.(id.eq.10)))
go to 10
352 IF(((id.eq.2).or.(id.eq.4)).or.(id.eq.10)) j=j+1
355 IF(yp(k).lt.yp(k-1)) in=0
357 IF(xp(k).gt.xp(k-1)) il=0
359 IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
361 107
FORMAT(/,
"***poINT same as poINT before last***",/)
362 IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
363 1.and.k.ge.3)
go to 2000
365 IF(kprt.eq.1)
WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
369 IF(k.eq.nt)
WRITE(66,106)
370 IF(k.eq.nt)
go to 2000
373 IF(jy.gt.(ny-2))
go to 9
375 CALL fit(2,y(jy-1),y(jy),y(jy+1),y(jy+2),p(ix,jy-1),p(ix,jy),
376 1 p(ix,jy+1),p(ix,jy+2),yp(k),v,yq)
378 CALL cubic(y(jy-1),y(jy),y(jy+1),y(jy+2),ws,a,b,c,f)
379 gsq(k)=a*grsq(ix,jy-1)+b*grsq(ix,jy)+c*grsq(ix,jy+1)
383 CALL fit(2,y(1),y(2),y(3),y(4),p(ix,1),p(ix,2),p(ix,3),p(ix,4),
386 CALL cubic(y(1),y(2),y(3),y(4),ws,a,b,c,f)
387 gsq(k)=a*grsq(ix,1)+b*grsq(ix,2)+c*grsq(ix,3)+f*grsq(ix,4)
390 CALL fit(2,y(ny-3),y(ny-2),y(ny-1),y(ny),p(ix,ny-3),p(ix,ny-2),
391 1 p(ix,ny-1),p(ix,ny),yp(k),v,yq)
393 CALL cubic(y(ny-3),y(ny-2),y(ny-1),y(ny),ws,a,b,c,f)
394 gsq(k)=a*grsq(ix,ny-3)+b*grsq(ix,ny-2)+c*grsq(ix,ny-1)
396 14
IF(yp(k).le.y(jy)) yp(k)=y(jy)+ff
397 IF(yp(k).ge.y(jy+1)) yp(k)=y(jy+1)-ff
399 IF((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
400 1.and.k.ne.1)
go to 10
401 19
IF(ix.eq.1.and.(((id.eq.5).or.(id.eq.9)).or.(id.eq.12)))
go to 10
402 IF(ix.eq.nx.and.(((id.eq.3).or.(id.eq.7)).or.(id.eq.11)))
go to 10
404 IF((id.eq.7).or.(id.eq.9)) j=j-1
406 IF(yp(k).lt.yp(k-1)) in=0
408 IF(xp(k).gt.xp(k-1)) il=0
409 IF(((id.eq.3).or.(id.eq.7.and.il.eq.0)).or.(id.eq.11.and.il.eq.0))
411 IF(kprt.eq.1)
WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
414 IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
416 IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
417 1.and.k.ge.3)
go to 2000
421 IF(k.eq.nt)
WRITE(66,106)
422 IF(k.eq.nt)
go to 2000
426 IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
428 IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
429 1.and.k.ge.3)
go to 2000
434 IF((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
435 1.and.k.ne.1)
go to 10
436 18
IF(((i.eq.1.or.i.eq.nx).and.(j.eq.1.or.j.eq.ny)).and.k.ne.1)
439 IF(yp(k).lt.yp(k-1)) in=0
441 IF(xp(k).gt.xp(k-1)) il=0
442 IF(kprt.eq.1)
WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
454 IF(abs(xp(j)-x1).gt.e)
go to 2008
455 IF(abs(yp(j)-y1).le.e)
go to 2030
469 WRITE(6,*)
'error in furpl--fur1'
471 WRITE(66,2002) v,(ixx,x(ixx),ixx=1,nx),(iyy,y(iyy),iyy=1,ny)
472 2002
FORMAT(
" furplm diagnostic center",/,
473 .
" value =",e14.6,/,
" x and y grids follow",/,
474 .(i5,e14.6,i5,e14.6,i5,e14.6,i5,e14.6))
476 2003
FORMAT(//,
" k ",
" xp yp p(i,j) ",
477 .
" i j ix jy ig in id iq il")
480 END SUBROUTINE furplm
481 SUBROUTINE cubic(p1,p2,p3,p4,p,a,b,c,d)
484 REAL(rprec) :: p1, p2, p3, p4, p, a, b, c, d, pp1, pp2, pp3, pp4
485 REAL(rprec) :: zero, p1p2, p2p3, p1p3, p1p4, p2p4, p3p4
486 INTEGER :: k12, k23, k34
490 k12=sign(1.5_dbl,p1p2)
491 IF(abs(p1p2).lt.zero)
go to 1
493 k23=sign(1.5_dbl,p2p3)
494 IF(abs(p2p3).lt.zero)
go to 2
496 IF(abs(p1p3).lt.zero)
go to 1
498 IF (abs(p1p4).lt.zero)
go to 3
500 IF (abs(p2p4).lt.zero)
go to 4
502 k34=sign(1.5_dbl,p3p4)
503 IF(abs(p3p4).lt.zero)
go to 4
504 IF (iabs(k12+k23+k34).eq.3)
go to 5
505 IF (iabs(k23+k34).eq.2)
go to 1
506 IF (iabs(k12+k23).eq.2)
go to 4
511 a=pp2*pp3*pp4/(p1p2*p1p3*p1p4)
512 b=-pp1*pp3*pp4/(p1p2*p2p3*p2p4)
513 c=pp1*pp2*pp4/(p1p3*p2p3*p3p4)
514 d=-pp1*pp2*pp3/(p1p4*p2p4*p3p4)
517 CALL parab(p2,p3,p4,p,b,c,d)
524 3
IF (iabs(k12+k23).eq.2)
go to 4
527 CALL parab(p1,p2,p3,p,a,b,c)
530 SUBROUTINE parab(p1,p2,p3,p,a,b,c)
533 REAL(rprec) :: p1, p2, p3, p, a, b, c, zero, yp, y, xp, x,
534 . p1p2, p2p3, p1p3, p1p4, p2p4, p3p4, pp1, pp2, pp3, pp4
535 INTEGER :: k12, k23, k34
544 k12=sign(1.5_dbl,p1p2)
545 IF(abs(p1p2).lt.zero)
go to 1
546 k23=sign(1.5_dbl,p2p3)
547 IF(abs(p2p3).lt.zero)
go to 2
548 IF(abs(p1p3).lt.zero)
go to 1
549 IF (iabs(k12+k23).ne.2)
go to 1
550 a=pp2*pp3/(p1p2*p1p3)
551 b=-pp1*pp3/(p1p2*p2p3)
552 c=pp1*pp2/(p1p3*p2p3)
563 SUBROUTINE fit(k,x1,x2,x3,x4,y1,y2,y3,y4,x,y,yp)
570 REAL(rprec) :: x1, x2, x3, x4, y1, y2, y3, y4, x, y, yp, f,
571 . c1, c2, c3, c4, d1, d2, d3, d4, d12, d13, d14, d23, d24, d34,
572 . crit, xa, xb, ya, yb, dydx
573 INTEGER :: k, iturn, i
575 c1=y1/((x1-x2)*(x1-x3)*(x1-x4))
576 c2=y2/((x2-x1)*(x2-x3)*(x2-x4))
577 c3=y3/((x3-x1)*(x3-x2)*(x3-x4))
578 c4=y4/((x4-x1)*(x4-x2)*(x4-x3))
590 f=(c1*d23+c2*d13+c3*d12)*d4+c4*d12*d3
591 yp=c1*(d23+d24+d34)+c2*(d13+d14+d34)
592 + +c3*(d12+d14+d24)+c4*(d12+d13+d23)
596 2
IF(y.ge.min(y1,y2,y3,y4).and.y.le.max(y1,y2,y3,y4))
GOTO 4
598 WRITE(6,11)x1,x2,x3,x4,y1,y2,y3,y4,y
601 500
FORMAT(
"error in sub fit at label 21")
602 WRITE(6,*)
'error in furpl--fit1'
605 11
FORMAT(
" fit",9e14.5)
606 4 crit=(abs(y1)+abs(y2)+abs(y3)+abs(y4))*1.e-05
612 IF((x-xa)*(x-xb).lt.0.)
go to 10
617 IF((x-xa)*(x-xb).lt.0.)
go to 10
622 IF((x-xa)*(x-xb).lt.0.)
go to 10
628 IF((y-ya)*(y-yb).lt.0.)
go to 1
634 IF((y-ya)*(y-yb).lt.0.)
go to 1
640 IF((y-ya)*(y-yb).lt.0.)
go to 1
644 13
IF(y.ne.y2)
go to 14
647 14
IF(y.ne.y3)
go to 15
650 15
IF(y.ne.y4)
go to 16
656 3
IF(abs(f-y).lt.crit) iturn=1
659 IF(abs(yp-dydx).lt..2*abs(yp))
go to 7
660 IF((f-y)*(ya-y).lt.0.)
go to 5
669 7
IF((f-y)*(ya-y).lt.0.)
go to 8
675 9 dydx=(yb-ya)/(xb-xa)
676 IF(abs(yp-dydx).lt..2*abs(yp))dydx=yp
678 IF(iturn.eq.1)
RETURN
681 10
IF((y-ya)*(y-yb).lt.0.)
go to 1