2 MODULE physical_constants
8 REAL(rprec),
PARAMETER :: pi=3.14159265358979323846264338328_rprec
9 REAL(rprec),
PARAMETER :: pio2=1.570796326794896619231321691_rprec
10 REAL(rprec),
PARAMETER :: twopi=6.28318530717958647692528677_rprec
11 REAL(rprec),
PARAMETER :: sqrt2=1.41421356237309504880168872_rprec
12 REAL(rprec),
PARAMETER :: degree = twopi / 360
13 REAL(rprec),
PARAMETER :: one = 1.0_rprec
14 REAL(rprec),
PARAMETER :: zero = 0.0_rprec
15 REAL(rprec),
PARAMETER :: two = 2.0_rprec
16 REAL(rprec),
PARAMETER :: half = 0.5_rprec
17 REAL(rprec),
PARAMETER :: three = 3.0_rprec
18 REAL(rprec),
PARAMETER :: four = 4.0_rprec
19 REAL(rprec),
PARAMETER :: six = 6.0_rprec
20 REAL(rprec),
PARAMETER :: twofour = 24.0_rprec
24 REAL(rprec),
PARAMETER :: mu0 = 2 * twopi * 1.0e-7_rprec
25 END MODULE physical_constants
26 subroutine furplm(p,grsq,x,y,v,xp,yp,gsq,np,nx,ny,ntx,nw,nh,is,js)
49 USE physical_constants
52 INTEGER :: js, is, ntx, ny, nx, np, nw, nh, nt, jst, nx2, ny2
53 INTEGER :: kprt, j, i, k, iq, ig, in, id, il, jy, ix, igp, l
54 INTEGER :: ixx, iyy, jdec
55 REAL(rprec) :: v, a, b, c, f, yq, dx, dy , dp, e, d, ff, ws
56 REAL(rprec) :: xp(*),yp(*), x1, y1
57 REAL(rprec) :: grsq(nw,nh),gsq(*)
58 REAL(rprec) :: p(nw,nh),x(nw),y(nh)
67 if(np.ne.-1)
go to 1999
71 1999
if(nt.gt.0)
go to 2004
82 e=10._dbl**(-1*4)*dx*three
94 13
if((v-p(i,j))*(p(i+1,j)-v).gt.zero)
go to 30
95 if(abs(p(i,j)-v).ge.d)
go to 32
103 if(jdec.lt.0)
go to 2005
107 2005
if(j.lt.1)
go to 10
110 30
if(i.gt.(nx-2))
go to 74
113 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),
114 1 p(i+2,j),xp(1),v,yq)
116 call cubic(x(i-1),x(i),x(i+1),x(i+2),ws,a,b,c,f)
117 gsq(1)=a*grsq(i-1,j)+b*grsq(i,j)+c*grsq(i+1,j)+f*grsq(i+2,j)
120 call fit(2,x(1),x(2),x(3),x(4),p(1,j),p(2,j),p(3,j),p(4,j),
123 call cubic(x(1),x(2),x(3),x(4),ws,a,b,c,f)
124 gsq(1)=a*grsq(1,j)+b*grsq(2,j)+c*grsq(3,j)+f*grsq(4,j)
127 call fit(2,x(nx-3),x(nx-2),x(nx-1),x(nx),p(nx-3,j),p(nx-2,j),
128 1 p(nx-1,j),p(nx,j),xp(1),v,yq)
130 call cubic(x(nx-3),x(nx-2),x(nx-1),x(nx),ws,a,b,c,f)
131 gsq(1)=a*grsq(nx-3,j)+b*grsq(nx-2,j)+c*grsq(nx-1,j)+f*grsq(nx,j)
134 if(xp(1).le.x(i)) xp(1)=x(i)+ff
135 if(xp(1).ge.x(i+1)) xp(1)=x(i+1)-ff
140 if(kprt.eq.1)
write(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
142 102
format(i5,3e13.6,10i5)
151 4
if(il.eq.1)
go to 21
155 if(i+1.gt.nx)
go to 10
159 if(abs(v-p(i+1,j+1)).ge.d)
go to 45
166 45
if(abs(v-p(i+1,j)).ge.d)
go to 33
175 33
if((v-p(i,j))*(p(i+1,j)-v).le.zero)
go to 50
180 50
if((v-p(i,j+1))*(p(i+1,j+1)-v).le.zero)
go to 51
185 51
if((v-p(i+1,j))*(p(i+1,j+1)-v).le.zero)
go to 52
194 if(i-1.lt.0)
go to 10
198 21
if(abs(v-p(i-1,j)).ge.d)
go to 46
204 46
if(abs(v-p(i-1,j+1)).ge.d)
go to 34
214 34
if((v-p(i,j+1))*(p(i-1,j+1)-v).le. zero)
go to 53
219 53
if((v-p(i-1,j+1))*(p(i-1,j)-v).le. zero)
go to 54
224 54
if((v-p(i,j))*(p(i-1,j)-v).le. zero)
go to 55
236 5
if(id.eq.15.and.il.eq.1) i=i-1
238 62
if(in.eq.1)
go to 12
242 if(j-1.lt.1)
go to 10
243 11
if(abs(v-p(i+1,j-1)).ge.d)
go to 42
250 42
if(abs(v-p(i,j-1)).ge.d)
go to 17
259 17
if((v-p(i+1,j-1))*(p(i+1,j)-v).le.zero)
go to 56
260 if(abs(v-p(i+1,j)).le.d)
go to 56
261 if(id.eq.15.and.igp.eq.7)
go to 56
266 56
if((v-p(i,j-1))*(p(i+1,j-1)-v).le.zero)
go to 57
271 57
if((v-p(i,j))*(p(i,j-1)-v).le.zero)
go to 44
272 if(abs(v-p(i,j)).le.d)
go to 44
277 44
if(iq.ne.1)
go to 48
288 12
if(j+1.gt.ny)
go to 10
293 if(abs(v-p(i+1,j+1)).ge.d)
go to 40
300 40
if(abs(v-p(i,j+1)).ge.d)
go to 61
309 61
if((v-p(i,j+1))*(p(i+1,j+1)-v).le. zero)
go to 59
314 59
if((v-p(i+1,j))*(p(i+1,j+1)-v).le. zero)
go to 60
315 if(abs(v-p(i+1,j)).le.d)
go to 60
320 60
if((v-p(i,j+1))*(p(i,j)-v).le. zero)
go to 41
321 if(abs(v-p(i,j)).le.d)
go to 41
326 41
if(iq.ne.1)
go to 49
340 if(k.eq.nt)
write(66,106)
341 106
format(/,
"***nt points found in furplm***",/)
342 if(k.eq.nt)
go to 2000
344 if(ix.gt.(nx-2))
go to 3
346 call fit(2,x(ix-1),x(ix),x(ix+1),x(ix+2),p(ix-1,jy),p(ix,jy),
347 1 p(ix+1,jy),p(ix+2,jy),xp(k),v,yq)
349 call cubic(x(ix-1),x(ix),x(ix+1),x(ix+2),ws,a,b,c,f)
351 gsq(k)=a*grsq(ix-1,jy)+b*grsq(ix,jy)+c*grsq(ix+1,jy)
356 call fit(2,x(1),x(2),x(3),x(4),p(1,jy),p(2,jy),p(3,jy),p(4,jy),
359 call cubic(x(1),x(2),x(3),x(4),ws,a,b,c,f)
360 gsq(k)=a*grsq(1,jy)+b*grsq(2,jy)+c*grsq(3,jy)+f*grsq(4,jy)
363 call fit(2,x(nx-3),x(nx-2),x(nx-1),x(nx),p(nx-3,jy),p(nx-2,jy),
364 1 p(nx-1,jy),p(nx,jy),xp(k),v,yq)
366 call cubic(x(nx-3),x(nx-2),x(nx-1),x(nx),ws,a,b,c,f)
367 gsq(k)=a*grsq(nx-3,jy)+b*grsq(nx-2,jy)+c*grsq(nx-1,jy)
370 if((id.eq.4).or.(id.eq.6)) i=i-1
371 if(xp(k).le.x(ix)) xp(k)=x(ix)+ff
372 if(xp(k).ge.x(ix+1)) xp(k)=x(ix+1)-ff
374 if((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
375 1.and.k.ne.1)
go to 10
376 20
if(jy.eq.1.and.(((id.eq.1).or.(id.eq.6)).or.(id.eq.8)))
go to 10
377 if(jy.eq.ny.and.(((id.eq.2).or.(id.eq.4)).or.(id.eq.10)))
go to 10
378 if(((id.eq.2).or.(id.eq.4)).or.(id.eq.10)) j=j+1
381 if(yp(k).lt.yp(k-1)) in=0
383 if(xp(k).gt.xp(k-1)) il=0
387 if((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
389 107
format(/,
"***point same as point before last***",/)
390 if((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
391 1.and.k.ge.3)
go to 2000
392 if(kprt.eq.1)
write(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
396 if(k.eq.nt)
write(66,106)
397 if(k.eq.nt)
go to 2000
400 if(jy.gt.(ny-2))
go to 9
402 call fit(2,y(jy-1),y(jy),y(jy+1),y(jy+2),p(ix,jy-1),p(ix,jy),
403 1 p(ix,jy+1),p(ix,jy+2),yp(k),v,yq)
405 call cubic(y(jy-1),y(jy),y(jy+1),y(jy+2),ws,a,b,c,f)
406 gsq(k)=a*grsq(ix,jy-1)+b*grsq(ix,jy)+c*grsq(ix,jy+1)
410 call fit(2,y(1),y(2),y(3),y(4),p(ix,1),p(ix,2),p(ix,3),p(ix,4),
413 call cubic(y(1),y(2),y(3),y(4),ws,a,b,c,f)
414 gsq(k)=a*grsq(ix,1)+b*grsq(ix,2)+c*grsq(ix,3)+f*grsq(ix,4)
417 call fit(2,y(ny-3),y(ny-2),y(ny-1),y(ny),p(ix,ny-3),p(ix,ny-2),
418 1 p(ix,ny-1),p(ix,ny),yp(k),v,yq)
420 call cubic(y(ny-3),y(ny-2),y(ny-1),y(ny),ws,a,b,c,f)
421 gsq(k)=a*grsq(ix,ny-3)+b*grsq(ix,ny-2)+c*grsq(ix,ny-1)
423 14
if(yp(k).le.y(jy)) yp(k)=y(jy)+ff
424 if(yp(k).ge.y(jy+1)) yp(k)=y(jy+1)-ff
426 if((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
427 1.and.k.ne.1)
go to 10
428 19
if(ix.eq.1.and.(((id.eq.5).or.(id.eq.9)).or.(id.eq.12)))
go to 10
429 if(ix.eq.nx.and.(((id.eq.3).or.(id.eq.7)).or.(id.eq.11)))
go to 10
431 if((id.eq.7).or.(id.eq.9)) j=j-1
433 if(yp(k).lt.yp(k-1)) in=0
435 if(xp(k).gt.xp(k-1)) il=0
436 if(((id.eq.3).or.(id.eq.7.and.il.eq.0)).or.(id.eq.11.and.il.eq.0))
438 if(kprt.eq.1)
write(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
443 if((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
445 if((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
446 1.and.k.ge.3)
go to 2000
449 if(k.eq.nt)
write(66,106)
450 if(k.eq.nt)
go to 2000
453 if((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
455 if((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
456 1.and.k.ge.3)
go to 2000
460 if((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
461 1.and.k.ne.1)
go to 10
462 18
if(((i.eq.1.or.i.eq.nx).and.(j.eq.1.or.j.eq.ny)).and.k.ne.1)
465 if(yp(k).lt.yp(k-1)) in=0
467 if(xp(k).gt.xp(k-1)) il=0
468 if(kprt.eq.1)
write(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
480 if(abs(xp(j)-x1).gt.e)
go to 2008
481 if(abs(yp(j)-y1).le.e)
go to 2030
495 write(6,*)
'error in furpl--fur1'
497 write(66,2002) v,(ixx,x(ixx),ixx=1,nx),(iyy,y(iyy),iyy=1,ny)
498 2002
format(
" furplm diagnostic center",/,
499 .
" value =",e14.6,/,
" x and y grids follow",/,
500 .(i5,e14.6,i5,e14.6,i5,e14.6,i5,e14.6))
502 2003
format(//,
" k ",
" xp yp p(i,j) ",
503 .
" i j ix jy ig in id iq il")
507 subroutine fit(k,x1,x2,x3,x4,y1,y2,y3,y4,x,y,yp)
513 USE physical_constants
516 INTEGER :: k, iturn, i
517 REAL(rprec) :: x1, x2, x3, x4, y1, y2, y3, y4, x, y, yp
518 REAL(rprec) :: c1, c2, c3, c4, d1, d2, d3, d4, crit, xa, ya
519 REAL(rprec) :: d12, d13, d14, d23, d24, d34, f, xb, yb, dydx
521 c1=y1/((x1-x2)*(x1-x3)*(x1-x4))
522 c2=y2/((x2-x1)*(x2-x3)*(x2-x4))
523 c3=y3/((x3-x1)*(x3-x2)*(x3-x4))
524 c4=y4/((x4-x1)*(x4-x2)*(x4-x3))
536 f=(c1*d23+c2*d13+c3*d12)*d4+c4*d12*d3
537 yp=c1*(d23+d24+d34)+c2*(d13+d14+d34)
538 + +c3*(d12+d14+d24)+c4*(d12+d13+d23)
542 2
if(y.ge.min(y1,y2,y3,y4).and.y.le.max(y1,y2,y3,y4))
go to 4
544 write(6,11)x1,x2,x3,x4,y1,y2,y3,y4,y
547 500
format(
"error in sub fit at label 21")
548 write(6,*)
'error in furpl--fit1'
551 11
format(
" fit",9e14.5)
552 4 crit=(abs(y1)+abs(y2)+abs(y3)+abs(y4))*1.d-05
558 if((x-xa)*(x-xb).lt.zero)
go to 10
563 if((x-xa)*(x-xb).lt.zero)
go to 10
568 if((x-xa)*(x-xb).lt.zero)
go to 10
574 if((y-ya)*(y-yb).lt.zero)
go to 1
580 if((y-ya)*(y-yb).lt.zero)
go to 1
586 if((y-ya)*(y-yb).lt.zero)
go to 1
590 13
if(y.ne.y2)
go to 14
593 14
if(y.ne.y3)
go to 15
596 15
if(y.ne.y4)
go to 16
602 3
if(abs(f-y).lt.crit) iturn=1
605 if(abs(yp-dydx).lt..2*abs(yp))
go to 7
606 if((f-y)*(ya-y).lt.zero)
go to 5
615 7
if((f-y)*(ya-y).lt.zero)
go to 8
621 9 dydx=(yb-ya)/(xb-xa)
622 if(abs(yp-dydx).lt..2*abs(yp))dydx=yp
624 if(iturn.eq.1)
return
627 10
if((y-ya)*(y-yb).lt.zero)
go to 1
630 subroutine cubic(p1,p2,p3,p4,p,a,b,c,d)
632 USE physical_constants
634 INTEGER :: k12, k23, k34
635 REAL(rprec) :: d, c, b, a, p4, p3, p2, p1
636 REAL(rprec) :: p1p2, p2p3, p1p3, p1p4, p2p4, p3p4
637 REAL(rprec) :: pp1, pp2, pp3, pp4, p
640 k12=sign(1.5_dbl,p1p2)
641 if(abs(p1p2).lt.zero)
go to 1
643 k23=sign(1.5_dbl,p2p3)
644 if(abs(p2p3).lt.zero)
go to 2
646 if(abs(p1p3).lt.zero)
go to 1
648 if (abs(p1p4).lt.zero)
go to 3
650 if (abs(p2p4).lt.zero)
go to 4
652 k34=sign(1.5_dbl,p3p4)
653 if(abs(p3p4).lt.zero)
go to 4
654 if (iabs(k12+k23+k34).eq.3)
go to 5
655 if (iabs(k23+k34).eq.2)
go to 1
656 if (iabs(k12+k23).eq.2)
go to 4
661 a=pp2*pp3*pp4/(p1p2*p1p3*p1p4)
662 b=-pp1*pp3*pp4/(p1p2*p2p3*p2p4)
663 c=pp1*pp2*pp4/(p1p3*p2p3*p3p4)
664 d=-pp1*pp2*pp3/(p1p4*p2p4*p3p4)
667 call parab(p2,p3,p4,p,b,c,d)
674 3
if (iabs(k12+k23).eq.2)
go to 4
677 call parab(p1,p2,p3,p,a,b,c)
680 subroutine parab(p1,p2,p3,p,a,b,c)
682 USE physical_constants
684 INTEGER :: k12, k23, k34
685 REAL(rprec) :: d, c, b, a, p4, p3, p2, p1
686 REAL(rprec) :: p1p2, p2p3, p1p3, p1p4, p2p4, p3p4
687 REAL(rprec) :: pp1, pp2, pp3, pp4, p
696 k12=sign(1.5_dbl,p1p2)
697 if(abs(p1p2).lt.zero)
go to 1
698 k23=sign(1.5_dbl,p2p3)
699 if(abs(p2p3).lt.zero)
go to 2
700 if(abs(p1p3).lt.zero)
go to 1
701 if (iabs(k12+k23).ne.2)
go to 1
702 a=pp2*pp3/(p1p2*p1p3)
703 b=-pp1*pp3/(p1p2*p2p3)
704 c=pp1*pp2/(p1p3*p2p3)