CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C subroutine subdipol C Now we want to subtract off the reference magnetic field "BNQF?". #include "param.inc" ****#include "meter.inc" #include "bzero.inc" #include "global_dims.inc" #include "mhd_P++_pointers.inc" #include "fortran_P++_arrays.inc" EXTERNAL BXqq0,BYqq0,BZqq0 #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in subdipol.F::subdipol(...)" #endif do k=1,nl_k do j=1,nl_j do i=1,nl_i+1 #ifdef BZEROS call g2int(bnqfi(i,j,k),bsqqi(i,j,k),bxqfi(i,j,k), $ byqfi(i,j,k),bzqfi(i,j,k),bxqni(i,j,k), $ byqni(i,j,k),bzqni(i,j,k),face_i(i,j,k), $ x(i,j,k),x(i,j+1,k),x(i,j,k+1),x(i,j+1,k+1), $ y(i,j,k),y(i,j+1,k),y(i,j,k+1),y(i,j+1,k+1), $ z(i,j,k),z(i,j+1,k),z(i,j,k+1),z(i,j+1,k+1)) #else bnqfi(i,j,k) = 0. bsqqi(i,j,k) = 0. bxqfi(i,j,k) = 0. byqfi(i,j,k) = 0. bzqfi(i,j,k) = 0. bxqni(i,j,k) = 0. byqni(i,j,k) = 0. bzqni(i,j,k) = 0. #endif enddo enddo enddo C do k=1,nl_k do j=1,nl_j+1 do i=1,nl_i #ifdef BZEROS call g2int(bnqfj(i,j,k),bsqqj(i,j,k),bxqfj(i,j,k), $ byqfj(i,j,k),bzqfj(i,j,k),bxqnj(i,j,k), $ byqnj(i,j,k),bzqnj(i,j,k),face_j(i,j,k), $ x(i,j,k),x(i,j,k+1),x(i+1,j,k),x(i+1,j,k+1), $ y(i,j,k),y(i,j,k+1),y(i+1,j,k),y(i+1,j,k+1), $ z(i,j,k),z(i,j,k+1),z(i+1,j,k),z(i+1,j,k+1)) if ( jlow_bound .and. j .eq. 1) bnqfj(i,j,k) = 0. if ( jhigh_bound .and. j .eq. nl_j+1 ) bnqfj(i,j,k) = 0. #else bnqfj(i,j,k) = 0. bsqqj(i,j,k) = 0. bxqfj(i,j,k) = 0. byqfj(i,j,k) = 0. bzqfj(i,j,k) = 0. bxqnj(i,j,k) = 0. byqnj(i,j,k) = 0. bzqnj(i,j,k) = 0. #endif enddo enddo enddo C do k=1,nl_k+1 do j=1,nl_j do i=1,nl_i #ifdef BZEROS call g2int(bnqfk(i,j,k),bsqqk(i,j,k),bxqfk(i,j,k), $ byqfk(i,j,k),bzqfk(i,j,k),bxqnk(i,j,k), $ byqnk(i,j,k),bzqnk(i,j,k),face_k(i,j,k), $ x(i,j,k),x(i+1,j,k),x(i,j+1,k),x(i+1,j+1,k), $ y(i,j,k),y(i+1,j,k),y(i,j+1,k),y(i+1,j+1,k), $ z(i,j,k),z(i+1,j,k),z(i,j+1,k),z(i+1,j+1,k)) #else bnqfk(i,j,k) = 0. bsqqk(i,j,k) = 0. bxqfk(i,j,k) = 0. byqfk(i,j,k) = 0. bzqfk(i,j,k) = 0. bxqnk(i,j,k) = 0. byqnk(i,j,k) = 0. bzqnk(i,j,k) = 0. #endif enddo enddo enddo C do 382 k=1,nl_k do 382 j=1,nl_j+1 do 382 i=1,nl_i bj_g(i,j,k) = bj_g(i,j,k) - bnqfj(i,j,k) 382 continue C do 384 k=1,nl_k do 384 j=1,nl_j do 384 i=1,nl_i+1 bi_g(i,j,k) = bi_g(i,j,k) - bnqfi(i,j,k) 384 continue C do 386 k=1,nl_k+1 do 386 j=1,nl_j do 386 i=1,nl_i bk_g(i,j,k) = bk_g(i,j,k) - bnqfk(i,j,k) 386 continue c C CALL BZZ(BX_g,BY_g,BZ_g,BI_g,BJ_g,BK_g) C DO 400 K=1,NL_K DO 400 J=1,NL_J+1 DO 400 I=1,NL_I+1 #ifdef BZEROS BZqK(I,J,K)=GINT(BZqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I,J,K+1),Y(I,J,K+1),Z(I,J,K+1)) BYqK(I,J,K)=GINT(BYqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I,J,K+1),Y(I,J,K+1),Z(I,J,K+1)) 400 BXqK(I,J,K)=GINT(BXqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I,J,K+1),Y(I,J,K+1),Z(I,J,K+1)) #else BZqK(I,J,K)=0. BYqK(I,J,K)=0. 400 BXqK(I,J,K)=0. #endif c c DO 410 K=1,NL_K+1 DO 410 J=1,NL_J DO 410 I=1,NL_I+1 #ifdef BZEROS BZqJ(I,J,K)=GINT(BZqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I,J+1,K),Y(I,J+1,K),Z(I,J+1,K)) BYqJ(I,J,K)=GINT(BYqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I,J+1,K),Y(I,J+1,K),Z(I,J+1,K)) 410 BXqJ(I,J,K)=GINT(BXqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I,J+1,K),Y(I,J+1,K),Z(I,J+1,K)) #else BZqJ(I,J,K)=0. BYqJ(I,J,K)=0. 410 BXqJ(I,J,K)=0. #endif DO 420 K=1,NL_K+1 DO 420 J=1,NL_J+1 DO 420 I=1,NL_I #ifdef BZEROS BZqI(I,J,K)=GINT(BZqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I+1,J,K),Y(I+1,J,K),Z(I+1,J,K)) BYqI(I,J,K)=GINT(BYqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I+1,J,K),Y(I+1,J,K),Z(I+1,J,K)) 420 BXqI(I,J,K)=GINT(BXqq0,X(I,J,K),Y(I,J,K), $ Z(I,J,K),X(I+1,J,K),Y(I+1,J,K),Z(I+1,J,K)) #else BZqI(I,J,K)=0. BYqI(I,J,K)=0. 420 BXqI(I,J,K)=0. #endif C return end C c subroutine mijset(i,j) #include "mdims.inc" #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in subdipol.F::mijset(...)" #endif mi =i mi2 = mi/2 mi2p1 = mi2+1 mip1 = mi+1 mip2 = mip1+1 mip3 = mip2+1 mip4 = mip3+1 mip5 = mip4+1 mim1 = mi-1 mim2 = mim1-1 mim3 = mim2-1 mim4 = mim3-1 C mj =j mj2 = mj/2 mj2p1 = mj2+1 mj2p2 = mj2+2 mjp1 = mj+1 mjp2 = mjp1+1 mjp3 = mjp2+1 mjp4 = mjp3+1 mjp5 = mjp4+1 mjm1 = mj-1 mjm2 = mjm1-1 mjm3 = mjm2-1 mjm4 = mjm3-1 C return end FUNCTION GINT(FX,X,Y,Z,X1,Y1,Z1) SAVE C Integrate FX(x,y,z) over the line (X,Y,Z) to (X1,Y1,Z1). C This subroutine does Gaussian integration with the first twelve C Legendre polynomials as the basis fuctions. C Abromowitz and Stegun page 916. C EXTERNAL FX DIMENSION WT(6),A(6) C The A() are the positive zeros of 12th order Legendre polynomial. DATA A/0.1252334085, 0.3678314989, 0.5873179542, $ 0.7699026741, 0.9041172563, 0.9815606342/ C Next are the Gaussian Integration coefficients for a C 12th order polynomial. DATA WT/0.2491470458, 0.2334925365, 0.2031674267, $ 0.1600783285, 0.1069393259, 0.0471753363/ C #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in subdipol.F::GINT(...)" #endif DX = (X1-X)*0.5 DY = (Y1-Y)*0.5 DZ = (Z1-Z)*0.5 XBAR = 0.5*(X+X1) YBAR = 0.5*(Y+Y1) ZBAR = 0.5*(Z+Z1) C SUM = 0. C DO 200 K=1,6 200 SUM = SUM + WT(K)*( $ FX(XBAR+A(K)*DX,YBAR+A(K)*DY,ZBAR+A(K)*DZ) + $ FX(XBAR-A(K)*DX,YBAR-A(K)*DY,ZBAR-A(K)*DZ) ) C GINT = 0.5*SUM C RETURN END C function fluxvol(r,L,dum) #include "param.inc" ! to get Rion real L C #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in subdipol.F::fluxvol(...)" #endif q = 1./(Rion*L) x = r*q fluxvol = sqrt((1.-3.*x/4.)/(1.-x))*r**3 C return end C C subroutine g2int(fnorm,FPRESS,fx,fy,fz,fnx,fny,fnz,face, $ x0,x1,x2,x3,y0,y1,y2,y3,z0,z1,z2,z3) SAVE C This subroutine does Gaussian integration of the function C C The Gaussian integration uses the first twelve Legendre polynomials C as the basis functions. Abromowitz and Stegun page 916. #ifndef T3E double precision bx,by,bz,bsq #endif real dx(3),dy(3),dz(3) real WT(12),A(12) DATA A/0.1252334085, 0.3678314989, 0.5873179542, $ 0.7699026741, 0.9041172563, 0.9815606342, $ -0.1252334085, -0.3678314989, -0.5873179542, $ -0.7699026741, -0.9041172563, -0.9815606342/ DATA WT/0.2491470458, 0.2334925365, 0.2031674267, $ 0.1600783285, 0.1069393259, 0.0471753363, $ 0.2491470458, 0.2334925365, 0.2031674267, $ 0.1600783285, 0.1069393259, 0.0471753363/ C #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in subdipol.F::g2int(...)" #endif dx(1) = x1-x0 dx(2) = x2-x0 dx(3) = x3+x0-x2-x1 dy(1) = y1-y0 dy(2) = y2-y0 dy(3) = y3+y0-y2-y1 dz(1) = z1-z0 dz(2) = z2-z0 dz(3) = z3+z0-z2-z1 C fnorm = 0. FPRESS = 0. fx = 0. fy = 0. fz = 0. fnx = 0. fny = 0. fnz = 0. C do 200 i=1,12 do 200 j=1,12 eta = 0.5*(1. + a(i)) psi = 0.5*(1. + a(j)) C x = x0 + dx(1)*eta + dx(2)*psi + dx(3)*eta*psi y = y0 + dy(1)*eta + dy(2)*psi + dy(3)*eta*psi z = z0 + dz(1)*eta + dz(2)*psi + dz(3)*eta*psi C xeta = dx(1) + dx(3)*psi yeta = dy(1) + dy(3)*psi zeta = dz(1) + dz(3)*psi C xpsi = dx(2) + dx(3)*eta ypsi = dy(2) + dy(3)*eta zpsi = dz(2) + dz(3)*eta C C note the fancy footwork with 1.e10 to make the vaxvms version behave C xn = 1.e-10*(yeta*zpsi-zeta*ypsi) yn = 1.e-10*(zeta*xpsi-xeta*zpsi) zn = 1.e-10*(xeta*ypsi-xpsi*yeta) bx = bxqq0(x,y,z) by = byqq0(x,y,z) bz = bzqq0(x,y,z) C bdotn = 2.5e9*(xn*bx + yn*by + zn*bz)*wt(i)*wt(j) rnsq = xn*xn+yn*yn+zn*zn rn = sqrt(rnsq) area = 2.5e9*rn*wt(i)*wt(j) BSQ = (BX*BX+BY*BY+BZ*BZ) fnorm = fnorm + bdotn fpress = fpress + area*BSQ fx = fx + area*bx fy = fy + area*by fz = fz + area*bz fnx = fnx + bdotn*bx fny = fny + bdotn*by fnz = fnz + bdotn*bz C 200 continue C fnx = fnx/(fnorm + 1.e-20) fny = fny/(fnorm + 1.e-20) fnz = fnz/(fnorm + 1.e-20) fnorm = fnorm/face fpress = fpress/face fx = fx/face fy = fy/face fz = fz/face C C return end C subroutine gauss(gpt,gwt) SAVE dimension gpt(12),gwt(12) DIMENSION WT(12),A(12) DATA A/0.1252334085, 0.3678314989, 0.5873179542, $ 0.7699026741, 0.9041172563, 0.9815606342, $ -0.1252334085, -0.3678314989, -0.5873179542, $ -0.7699026741, -0.9041172563, -0.9815606342/ DATA WT/0.2491470458, 0.2334925365, 0.2031674267, $ 0.1600783285, 0.1069393259, 0.0471753363, $ 0.2491470458, 0.2334925365, 0.2031674267, $ 0.1600783285, 0.1069393259, 0.0471753363/ #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in subdipol.F::gauss(...)" #endif do 200 i=1,6 gwt(i) = 0.5*wt(7-i) gpt(i) = 0.5 - 0.5*a(7-i) gwt(6+i) = 0.5*wt(i) gpt(6+i) = 0.5 + 0.5*a(i) 200 continue return end C