C---------------------------------------------------------------------- C Function performs red-black Gauss-Seidel relaxation. C C Date of last modification: C Tue Feb 20 14:40:47 PST 2001 C---------------------------------------------------------------------- double precision function res_relax( alpha, a, c, br, bz, & krr, phi, pi, phi_r, & phi_z, Nr, Nz, r, z, & order, mg_rhs1, mg_rhs2, & mg_rhs3, mg_rhs4) C---------------------------------------------------------------------- C Declare variables (explicitly). C---------------------------------------------------------------------- implicit none logical ltrace parameter ( ltrace = .false. ) logical ltrace0 parameter ( ltrace0 = .false. ) logical ltrace1 parameter ( ltrace1 = .false. ) character*10 cdnm parameter ( cdnm = 'res_relax' ) integer Nr, Nz integer i, j, k integer l, m integer isw, jsw, pass integer maxpass, step real*8 alpha(Nr, Nz) real*8 a(Nr,Nz), c(Nr, Nz) real*8 br(Nr,Nz), bz(Nr, Nz) real*8 krr(Nr, Nz) real*8 phi(Nr, Nz), pi(Nr, Nz) real*8 phi_r(Nr, Nz), phi_z(Nr, Nz) real*8 r(Nr), z(Nz) real*8 dr, dz real*8 s1, s2, s3, s4, s5, s6 real*8 Jac(4,4) real*8 res(4), lhs(4) real*8 du(4) real*8 cpi parameter ( cpi = 3.141 5926 5358 9793 d0 ) integer info integer W, B integer order real*8 mg_rhs1(Nr, Nz) real*8 mg_rhs2(Nr, Nz) real*8 mg_rhs3(Nr, Nz) real*8 mg_rhs4(Nr, Nz) C====================================================================== C====================================================================== C---------------------------------------------------------------------- C Some initializations and tracing. C---------------------------------------------------------------------- res_relax = 0.0d0 dr = r(2) - r(1) dz = z(2) - z(1) if ( ltrace ) then write(0,*) cdnm,': Nr =', Nr write(0,*) cdnm,': Nz =', Nz write(0,*) cdnm,': dr =', dr write(0,*) cdnm,': dz =', dz write(0,*) cdnm,': order =', order endif C---------------------------------------------------------------------- C Set the boundary conditions at r=0 for a, br, bzi. C---------------------------------------------------------------------- do j = 1, Nz a(1,j) = 1.0d0 br(1,j) = 0.0d0 bz(1,j) = 0.0d0 end do C---------------------------------------------------------------------- C Implement the red-black ralaxation. C---------------------------------------------------------------------- jsw = 1 if (order .eq. 1) then maxpass = 2 step = 2 else if (order .eq. 0 ) then maxpass = 1 step = 1 endif do pass = 1, maxpass isw = jsw do j = 1, Nz do i = isw + 1, Nr - 1, step C---------------------------------------------------------------------- C Calculation of the Jacobian. (Values/definitions of s1, ..., s6, C and the elements of array Jac are stored in file jacob.) C---------------------------------------------------------------------- C---------------------------------------------------------------------- C Set the integers W, B to take into account the periodic C boundary conditions. C---------------------------------------------------------------------- W = 1 B = 1 if ( j .eq. 1) B = j - Nz if ( j .eq. Nz) W = -j + 1 C---------------------------------------------------------------------- C Define s1, ..., s6, and the elements of arrady Jac from file C jacob. C---------------------------------------------------------------------- s2 = 1.D0/2.D0 s5 = 4.D0*a(i,j)*c(i,j)**3*dz**2*dr**2+16.D0*c(i,j)**3*r #(i)**2*cpi*a(i,j)*phi_r(i,j)**2*dz**2*dr**2+4.D0*c(i,j)*r #(i)**2*a(i,j)**2*a(i,j+W)*dr**2+4.D0*a(i,j)*r(i)*c(i,j) #**2*dr*dz**2*c(i+1,j)-4.D0*a(i,j)*r(i)*c(i,j)**2*dr*dz** #2*c(i-1,j)-r(i)**2*a(i,j)**2*dr**2*a(i,j+W)*c(i,j+W)+r #(i)**2*a(i,j)**2*dr**2*a(i,j+W)*c(i,j-B)-4.D0*c(i,j)**3. #D0*r(i)*dz**2*a(i+1,j)*dr+c(i,j)**2*r(i)**2*dz**2*a(i-1 #,j)*c(i+1,j)-c(i,j)**2*r(i)**2*dz**2*a(i-1,j)*c(i-1,j) s4 = s5+4.D0*c(i,j)**3*r(i)*dz**2*a(i-1,j)*dr+r(i)**2*a #(i,j)**2*dr**2*a(i,j-B)*c(i,j+W)-r(i)**2*a(i,j)**2*dr #**2*a(i,j-B)*c(i,j-B)+4.D0*c(i,j)*r(i)**2*a(i,j)**2*a( #i,j-B)*dr**2-c(i,j)**2*r(i)**2*dz**2*a(i+1,j)*c(i+1,j)+c #(i,j)**2*r(i)**2*dz**2*a(i+1,j)*c(i-1,j)-4.D0*c(i,j)**3 #.D0*r(i)**2*a(i,j)*dz**2-12.D0*a(i,j)**3*r(i)**2*c(i #,j)*dr**2-4.D0*a(i,j)**3*c(i,j)**3*dz**2*dr**2+16.D0*a(i, #j)**3*r(i)**2*cpi*c(i,j)*phi_z(i,j)**2*dz**2*dr**2 s5 = 1.D0/dr**2 s3 = s4*s5 s1 = s2*s3 s2 = 1.D0/c(i,j)**3/a(i,j)**3/r(i)**2/dz**2 Jac(1,1) = s1*s2 s2 = -1.D0/4.D0 s6 = c(i,j)*r(i)**2*a(i,j)**2*dr**2*a(i,j+W)*alpha(i,j+ #W)-2.D0*r(i)**2*a(i,j)**2*alpha(i,j)*dr**2*a(i,j-B)*c(i #,j-B)+2.D0*r(i)**2*a(i,j)**2*alpha(i,j)*dr**2*a(i,j-B)*c #(i,j+W)+8.D0*c(i,j)*r(i)**2*a(i,j)**2*alpha(i,j)*a(i,j #-B)*dr**2+8.D0*c(i,j)*r(i)**2*a(i,j)**2*alpha(i,j)*a(i, #j+W)*dr**2+16.D0*a(i,j)*r(i)*c(i,j)**2*alpha(i,j)*dr*dz**2* #c(i+1,j)-16.D0*a(i,j)*r(i)*c(i,j)**2*alpha(i,j)*dr*dz**2*c #(i-1,j)-16.D0*c(i,j)**3*r(i)**2*a(i,j)*dz**2*alpha(i,j) s5 = s6+2.D0*a(i,j)*r(i)**2*c(i,j)**2*dz**2*alpha(i-1,j) #*c(i-1,j)-2.D0*a(i,j)*r(i)**2*c(i,j)**2*dz**2*alpha(i-1 #,j)*c(i+1,j)-3.D0*c(i,j)**3*r(i)**2*dz**2*a(i+1,j)*alpha #(i+1,j)+3.D0*c(i,j)**3*r(i)**2*dz**2*a(i+1,j)*alpha(i-1 #,j)+2.D0*r(i)**2*a(i,j)**2*alpha(i,j)*dr**2*a(i,j+W)*c( #i,j-B)-2.D0*r(i)**2*a(i,j)**2*alpha(i,j)*dr**2*a(i,j+W)* #c(i,j+W)+24.D0*c(i,j)**3*r(i)*alpha(i,j)*dz**2*a(i-1,j)*dr #+6.D0*c(i,j)**2*r(i)**2*alpha(i,j)*dz**2*a(i-1,j)*c(i+1 #,j) s6 = s5-6.D0*c(i,j)**2*r(i)**2*alpha(i,j)*dz**2*a(i-1,j) #*c(i-1,j)-24.D0*c(i,j)**3*r(i)*alpha(i,j)*dz**2*a(i+1,j)*d #r-6.D0*c(i,j)**2*r(i)**2*alpha(i,j)*dz**2*a(i+1,j)*c(i+ #1,j)+6.D0*c(i,j)**2*r(i)**2*alpha(i,j)*dz**2*a(i+1,j)*c #(i-1,j)+8.D0*a(i,j)*c(i,j)**3*r(i)*dr*dz**2*alpha(i+1,j)+3. #D0*c(i,j)**3*r(i)**2*dz**2*a(i-1,j)*alpha(i+1,j)+2.D0*a #(i,j)*r(i)**2*c(i,j)**2*dz**2*alpha(i+1,j)*c(i+1,j)+c(i #,j)*r(i)**2*a(i,j)**2*dr**2*a(i,j-B)*alpha(i,j-B) s4 = s6-c(i,j)*r(i)**2*a(i,j)**2*dr**2*a(i,j+W)*alpha(i #,j-B)-8.D0*a(i,j)*c(i,j)**3*r(i)*dr*dz**2*alpha(i-1,j)+8.D0 #*c(i,j)**3*r(i)**2*a(i,j)*dz**2*alpha(i-1,j)+8.D0*c(i,j #)**3*r(i)**2*a(i,j)*dz**2*alpha(i+1,j)+16.D0*a(i,j)*c(i #,j)**3*alpha(i,j)*dz**2*dr**2-c(i,j)*r(i)**2*a(i,j)**2.D #0*dr**2*a(i,j-B)*alpha(i,j+W)+64.D0*c(i,j)**3*r(i)**2*al #pha(i,j)*cpi*a(i,j)*phi_r(i,j)**2*dz**2*dr**2-3.D0*c(i,j)**3 #.D0*r(i)**2*dz**2*a(i-1,j)*alpha(i-1,j)-2.D0*a(i,j)*r(i)** #2.D0*c(i,j)**2*dz**2*alpha(i+1,j)*c(i-1,j) s5 = 1.D0/dr**2 s3 = s4*s5 s1 = s2*s3 s2 = 1.D0/a(i,j)**4/c(i,j)**3/r(i)**2/dz**2 Jac(2,1) = s1*s2 Jac(3,1) = 0.D0 Jac(4,1) = 0.D0 Jac(1,2) = 1.D0/c(i,j)**2/a(i,j)**2*(r(i)**2*a(i,j)**4.D #0*dr**2*br(i,j-B)**2+r(i)**2*a(i,j)**4*dr**2*br(i,j+W #)**2+48.D0*br(i,j)**2*a(i,j)**2*c(i,j)**2*dz**2*dr* #*2+2.D0*r(i)**2*c(i,j)**2*a(i,j)**2*dz*br(i,j+W)*dr*b #z(i+1,j)+2.D0*r(i)**2*c(i,j)**2*a(i,j)**2*dz*br(i,j- #B)*dr*bz(i-1,j)-2.D0*r(i)**2*c(i,j)**2*a(i,j)**2*dz*b #r(i,j-B)*dr*bz(i+1,j)-2.D0*r(i)**2*c(i,j)**2*a(i,j)**2. #D0*dz*br(i,j+W)*dr*bz(i-1,j)-2.D0*r(i)**2*c(i,j)**4*dz** #2*bz(i+1,j)*bz(i-1,j)+r(i)**2*c(i,j)**4*dz**2*bz(i+1,j) #**2+r(i)**2*c(i,j)**4*dz**2*bz(i-1,j)**2-2.D0*r(i) #**2*a(i,j)**4*dr**2*br(i,j+W)*br(i,j-B)+16.D0*r(i)*krr( #i,j)*c(i,j)**2*br(i,j)*alpha(i,j)*dz**2*dr**2)/dr**2/dz**2/r #(i)**2/alpha(i,j)**3/4.D0 s4 = r(i)**2*c(i,j)**5*a(i,j)**2*dz**2*bz(i+1,j)**2.D #0/4.D0+r(i)**2*c(i,j)**5*a(i,j)**2*dz**2*bz(i-1,j)**2 #.D0/4.D0-r(i)**2*c(i,j)*a(i,j)**6*dr**2*br(i,j-B)**2/ #4.D0-3.D0/2.D0*c(i,j)**2*r(i)**2*alpha(i,j)**2*a(i,j) #*dz**2*a(i+1,j)*c(i+1,j)+3.D0/2.D0*c(i,j)**2*r(i)**2*alp #ha(i,j)**2*a(i,j)*dz**2*a(i+1,j)*c(i-1,j)-6.D0*c(i,j)**3.D #0*r(i)*alpha(i,j)**2*a(i,j)*dz**2*a(i+1,j)*dr s3 = s4+2.D0*a(i,j)**3*c(i,j)*r(i)**2*alpha(i,j)**2*a #(i,j-B)*dr**2+3.D0/2.D0*c(i,j)**2*r(i)**2*alpha(i,j)**2. #D0*a(i,j)*dz**2*a(i-1,j)*c(i+1,j)-3.D0/2.D0*c(i,j)**2*r(i) #**2*alpha(i,j)**2*a(i,j)*dz**2*a(i-1,j)*c(i-1,j)+8.D0*r #(i)**2*krr(i,j)**2*c(i,j)**3*alpha(i,j)**2*dz**2*dr #**2+6.D0*c(i,j)**3*r(i)*alpha(i,j)**2*a(i,j)*dz**2*a(i- #1,j)*dr+16.D0*cpi*a(i,j)**2*r(i)**2*c(i,j)**3*alpha(i #,j)**2*phi_r(i,j)**2*dz**2*dr**2 s4 = s3-a(i,j)**3*r(i)**2*alpha(i,j)**2*dr**2*a(i,j+W #)*c(i,j+W)/2.D0+a(i,j)**3*r(i)**2*alpha(i,j)**2*dr**2 #*a(i,j+W)*c(i,j-B)/2.D0+a(i,j)**3*r(i)**2*alpha(i,j)**2 #.D0*dr**2*a(i,j-B)*c(i,j+W)/2.D0-a(i,j)**3*r(i)**2*alpha #(i,j)**2*dr**2*a(i,j-B)*c(i,j-B)/2.D0+4.D0*a(i,j)**2*c( #i,j)**2*r(i)*dr*alpha(i,j)**2*dz**2*c(i+1,j) s2 = s4-4.D0*a(i,j)**2*c(i,j)**2*r(i)*dr*alpha(i,j)**2.D #0*dz**2*c(i-1,j)+r(i)**2*c(i,j)*a(i,j)**6*dr**2*br(i,j+ #W)*br(i,j-B)/2.D0+4.D0*a(i,j)**2*c(i,j)**3*alpha(i,j)**2 #.D0*dz**2*dr**2-r(i)**2*c(i,j)**5*a(i,j)**2*dz**2*bz( #i+1,j)*bz(i-1,j)/2.D0+2.D0*a(i,j)**3*c(i,j)*r(i)**2*alph #a(i,j)**2*a(i,j+W)*dr**2+8.D0*r(i)*krr(i,j)*c(i,j)**3*b #r(i,j)*a(i,j)**2*alpha(i,j)*dz**2*dr**2-r(i)**2*c(i,j)* #a(i,j)**6*dr**2*br(i,j+W)**2/4.D0 s3 = 1.D0/a(i,j)**5/r(i)**2 s1 = s2*s3 s2 = 1.D0/c(i,j)**3/alpha(i,j)**2/dz**2/dr**2 Jac(2,2) = s1*s2 Jac(3,2) = -4.D0/a(i,j)**2*(r(i)*krr(i,j)*alpha(i,j)+3.D0*br( #i,j)*a(i,j)**2)/r(i)**2/alpha(i,j)**2 Jac(4,2) = 0.D0 s4 = r(i)**2*c(i,j)*a(i,j)**3*alpha(i,j)*dr*br(i,j)+r( #i)**2*c(i,j)**3*a(i,j)*alpha(i,j)*dz*bz(i+1,j-B)/8.D0+c #(i,j)**2*alpha(i,j)*r(i)*a(i,j)*br(i,j)*dz**2*c(i-1,j)+2.D #0*c(i,j)**3*alpha(i,j)*a(i,j)*br(i,j)*dr*dz**2-c(i,j)**2.D #0*alpha(i,j)*r(i)*a(i,j)*br(i,j)*dz**2*c(i+1,j)+r(i)**2*c #(i,j)**3*a(i,j)*alpha(i,j)*dz*bz(i-1,j+W)/8.D0-r(i)**2* #c(i,j)**3*a(i,j)*alpha(i,j)*dz*bz(i+1,j+W)/8.D0-r(i)**2 #*c(i,j)**3*a(i,j)*alpha(i,j-B)*dz*bz(i+1,j)/4.D0 s3 = s4-r(i)**2*c(i,j)**3*alpha(i,j)*a(i,j+W)*dz*bz(i+1 #,j)/8.D0+r(i)**2*c(i,j)**3*alpha(i,j)*a(i,j-B)*dz*bz(i+ #1,j)/8.D0+r(i)**2*alpha(i,j)*dz*a(i,j)*c(i,j)**2*bz(i+1 #,j)*c(i,j-B)/8.D0-r(i)**2*alpha(i,j)*dz*a(i,j)*c(i,j)**2.D #0*bz(i+1,j)*c(i,j+W)/8.D0+r(i)**2*c(i,j)**3*a(i,j)*alph #a(i,j+W)*dz*bz(i+1,j)/4.D0-r(i)**2*c(i,j)**3*a(i,j)*alp #ha(i,j)*dz*bz(i-1,j-B)/8.D0-r(i)**2*c(i,j)**3*a(i,j)*al #pha(i,j+W)*dz*bz(i-1,j)/4.D0-r(i)**2*c(i,j)**3*alpha(i, #j)*a(i,j-B)*dz*bz(i-1,j)/8.D0-r(i)**2*alpha(i,j)*dz*a(i,j) #*c(i,j)**2*bz(i-1,j)*c(i,j-B)/8.D0 s4 = s3+r(i)**2*alpha(i,j)*dz*a(i,j)*c(i,j)**2*bz(i-1,j #)*c(i,j+W)/8.D0+r(i)**2*c(i,j)**3*a(i,j)*alpha(i,j-B)*d #z*bz(i-1,j)/4.D0+r(i)**2*c(i,j)**3*alpha(i,j)*a(i,j+W)* #dz*bz(i-1,j)/8.D0-r(i)**2*c(i,j)*a(i,j)**3*alpha(i,j)*d #r*br(i,j+W)/2.D0+r(i)**2*alpha(i,j)*a(i,j)**3*dr*br(i,j #-B)*c(i,j-B)/8.D0-r(i)**2*c(i,j)*a(i,j)**3*alpha(i,j)*d #r*br(i,j-B)/2.D0-3.D0/8.D0*r(i)**2*c(i,j)*alpha(i,j)*a(i,j #)**2*dr*a(i,j-B)*br(i,j-B)-r(i)**2*alpha(i,j)*a(i,j)**3 #.D0*dr*br(i,j-B)*c(i,j+W)/8.D0 s2 = s4-r(i)**2*c(i,j)*a(i,j)**3*dr*alpha(i,j+W)*br(i,j #-B)/4.D0+3.D0/8.D0*r(i)**2*c(i,j)*alpha(i,j)*a(i,j)**2*d #r*a(i,j+W)*br(i,j-B)+r(i)**2*c(i,j)*a(i,j)**3*dr*alpha #(i,j-B)*br(i,j-B)/4.D0-r(i)**2*alpha(i,j)*a(i,j)**3*dr*b #r(i,j+W)*c(i,j-B)/8.D0+r(i)**2*alpha(i,j)*a(i,j)**3*dr* #br(i,j+W)*c(i,j+W)/8.D0-3.D0/8.D0*r(i)**2*c(i,j)*alpha(i,j #)*a(i,j)**2*dr*a(i,j+W)*br(i,j+W)-r(i)**2*c(i,j)*a(i,j #)**3*dr*alpha(i,j-B)*br(i,j+W)/4.D0+r(i)**2*c(i,j)*a(i, #j)**3*dr*alpha(i,j+W)*br(i,j+W)/4.D0+3.D0/8.D0*r(i)**2*c #(i,j)*alpha(i,j)*a(i,j)**2*dr*a(i,j-B)*br(i,j+W) s3 = 1.D0/a(i,j)/dz**2 s1 = s2*s3 s2 = 1.D0/dr/r(i)**2/alpha(i,j)**3/c(i,j)**3 Jac(1,3) = s1*s2 s4 = -r(i)*c(i,j)*a(i,j)**5*dr*alpha(i,j-B)*br(i,j-B)/4.D0 #+r(i)*c(i,j)**3*alpha(i,j)*a(i,j)**2*a(i,j-B)*dz*bz(i+ #1,j)/8.D0+r(i)*c(i,j)*a(i,j)**5*dr*alpha(i,j+W)*br(i,j-B)/ #4.D0-3.D0/8.D0*r(i)*c(i,j)*alpha(i,j)*a(i,j)**4*dr*a(i,j-B #)*br(i,j+W)-3.D0/8.D0*r(i)*c(i,j)*alpha(i,j)*a(i,j)**4*dr* #a(i,j+W)*br(i,j-B)-r(i)*c(i,j)**3*alpha(i,j)*a(i,j)**2 #*a(i,j+W)*dz*bz(i+1,j)/8.D0 s3 = s4+3.D0/8.D0*r(i)*c(i,j)*alpha(i,j)*a(i,j)**4*dr*a(i, #j-B)*br(i,j-B)+r(i)*c(i,j)*a(i,j)**5*dr*alpha(i,j-B)*br(i #,j+W)/4.D0+r(i)*c(i,j)*a(i,j)**5*alpha(i,j)*dr*br(i,j+W)-r #(i)*c(i,j)**3*alpha(i,j)*a(i,j)**2*a(i,j-B)*dz*bz(i-1, #j)/8.D0-2.D0*r(i)*c(i,j)*a(i,j)**5*alpha(i,j)*dr*br(i,j)+3 #.D0/8.D0*r(i)*c(i,j)*alpha(i,j)*a(i,j)**4*dr*a(i,j+W)*br( #i,j+W) s4 = s3+r(i)*c(i,j)*a(i,j)**5*alpha(i,j)*dr*br(i,j-B)-r(i #)*c(i,j)*a(i,j)**5*dr*alpha(i,j+W)*br(i,j+W)/4.D0+r(i)*c( #i,j)**3*alpha(i,j)*a(i,j)**2*a(i,j+W)*dz*bz(i-1,j)/8.D0+ #r(i)*alpha(i,j)*a(i,j)**5*dr*br(i,j-B)*c(i,j+W)/4.D0-r(i) #*alpha(i,j)*a(i,j)**5*dr*br(i,j-B)*c(i,j-B)/4.D0 s2 = s4-r(i)*alpha(i,j)*a(i,j)**5*dr*br(i,j+W)*c(i,j+W)/4. #D0+r(i)*alpha(i,j)*a(i,j)**5*dr*br(i,j+W)*c(i,j-B)/4.D0-4. #D0*c(i,j)**3*krr(i,j)*alpha(i,j)**2*a(i,j)*dr*dz**2-3.D0 #*c(i,j)**3*r(i)*alpha(i,j)**2*krr(i,j)*a(i-1,j)*dz**2+3 #.D0*c(i,j)**3*r(i)*alpha(i,j)**2*krr(i,j)*a(i+1,j)*dz** #2-2.D0*r(i)*krr(i,j)*c(i,j)**2*alpha(i,j)**2*a(i,j)*dz* #*2*c(i+1,j)+2.D0*r(i)*krr(i,j)*c(i,j)**2*alpha(i,j)**2* #a(i,j)*dz**2*c(i-1,j) s3 = 1.D0/r(i)/dz**2 s1 = s2*s3 s2 = 1.D0/dr/alpha(i,j)**2/a(i,j)**4/c(i,j)**3 Jac(2,3) = s1*s2 Jac(3,3) = -(-c(i,j)*r(i)*dz**2*c(i+1,j)+c(i,j)*r(i)*dz**2*c(i- #1,j)+2.D0*c(i,j)**2*dz**2*dr+r(i)**2*a(i,j)**2*dr)/c( #i,j)**2/alpha(i,j)/dr/dz**2/r(i)**2 Jac(4,3) = 0.D0 s4 = -r(i)*dr*alpha(i,j)*a(i,j)**3*br(i,j+W)*c(i+1,j)/8.D0 #+r(i)*a(i,j)**3*c(i,j)*alpha(i+1,j)*br(i,j+W)*dr/4.D0+r(i #)*dr*alpha(i,j)*a(i,j)**3*br(i,j+W)*c(i-1,j)/8.D0+r(i)*alp #ha(i,j)*a(i,j)*c(i,j)**3*dz*bz(i,j)-r(i)*alpha(i,j)*a(i, #j)**3*c(i,j)*br(i-1,j-B)*dr/8.D0+r(i)*alpha(i,j)*a(i,j)**3 #.D0*c(i,j)*br(i+1,j-B)*dr/8.D0-r(i)*alpha(i,j)*a(i,j)**3*c #(i,j)*br(i+1,j+W)*dr/8.D0+r(i)*alpha(i,j)*a(i,j)**3*c(i,j #)*br(i-1,j+W)*dr/8.D0-r(i)*c(i,j)*alpha(i,j)*a(i,j)**2*a( #i+1,j)*br(i,j+W)*dr/8.D0 s3 = s4+r(i)*c(i,j)*alpha(i,j)*a(i,j)**2*a(i-1,j)*br(i,j+ #W)*dr/8.D0-r(i)*alpha(i,j)*a(i,j)*c(i,j)**3*dz*bz(i-1,j)/2 #.D0+r(i)*a(i,j)*c(i,j)**3*dz*alpha(i-1,j)*bz(i-1,j)/4.D0-r #(i)*a(i,j)**3*c(i,j)*alpha(i-1,j)*br(i,j+W)*dr/4.D0+alpha #(i,j)*a(i,j)*c(i,j)**3*dr*dz*bz(i-1,j)/2.D0-r(i)*alpha(i,j #)*c(i,j)**3*a(i+1,j)*dz*bz(i-1,j)/8.D0+r(i)*alpha(i,j)*c( #i,j)**3*a(i-1,j)*dz*bz(i-1,j)/8.D0-c(i,j)*alpha(i,j)*a(i,j #)**2*br(i,j)*dr**2*a(i,j-B)-alpha(i,j)*a(i,j)**3*c(i,j) #*dr**2*br(i,j-B)/2.D0+2.D0*a(i,j)**3*c(i,j)*br(i,j)*dr**2*a #lpha(i,j-B) s4 = s3-3.D0/8.D0*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i #-1,j)*c(i-1,j)+3.D0/8.D0*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2* #dz*bz(i-1,j)*c(i+1,j)+c(i,j)*alpha(i,j)*a(i,j)**2*br(i,j) #*dr**2*a(i,j+W)-2.D0*a(i,j)**3*c(i,j)*br(i,j)*dr**2*alpha( #i,j+W)+r(i)*a(i,j)*c(i,j)**3*dz*alpha(i+1,j)*bz(i+1,j)/4.D #0-r(i)*a(i,j)*c(i,j)**3*dz*alpha(i+1,j)*bz(i-1,j)/4.D0-r( #i)*alpha(i,j)*c(i,j)**3*a(i-1,j)*dz*bz(i+1,j)/8.D0-r(i)*al #pha(i,j)*a(i,j)*c(i,j)**3*dz*bz(i+1,j)/2.D0+3.D0/8.D0*r(i) #*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i+1,j)*c(i-1,j) s5 = s4+r(i)*alpha(i,j)*c(i,j)**3*a(i+1,j)*dz*bz(i+1,j)/8. #D0-3.D0/8.D0*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i+1,j) #*c(i+1,j)-r(i)*c(i,j)*alpha(i,j)*a(i,j)**2*a(i-1,j)*br(i #,j-B)*dr/8.D0+r(i)*c(i,j)*alpha(i,j)*a(i,j)**2*a(i+1,j)*br #(i,j-B)*dr/8.D0-r(i)*dr*alpha(i,j)*a(i,j)**3*br(i,j-B)*c( #i-1,j)/8.D0 s2 = s5+r(i)*dr*alpha(i,j)*a(i,j)**3*br(i,j-B)*c(i+1,j)/8. #D0+r(i)*a(i,j)**3*c(i,j)*alpha(i-1,j)*br(i,j-B)*dr/4.D0-r #(i)*a(i,j)**3*c(i,j)*alpha(i+1,j)*br(i,j-B)*dr/4.D0+alpha( #i,j)*a(i,j)**3*c(i,j)*dr**2*br(i,j+W)/2.D0-alpha(i,j)*a(i, #j)*c(i,j)**3*dr*dz*bz(i+1,j)/2.D0-r(i)*a(i,j)*c(i,j)**3 #*dz*alpha(i-1,j)*bz(i+1,j)/4.D0 s3 = 1.D0/dr**2/dz s1 = s2*s3 s2 = 1.D0/c(i,j)/r(i)/alpha(i,j)**3/a(i,j)**3 Jac(4,1) = s1*s2 s2 = -1.D0/8.D0 s6 = -r(i)*alpha(i,j)*a(i,j)**2*a(i-1,j)*br(i,j+W)*dr+8.D0 #*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i-1,j)+r(i)*alpha #(i,j)*a(i,j)**2*a(i+1,j)*br(i,j+W)*dr+3.D0*r(i)*alpha(i,j #)*c(i,j)**2*a(i+1,j)*dz*bz(i-1,j)+3.D0*r(i)*alpha(i,j)*c( #i,j)**2*a(i-1,j)*dz*bz(i+1,j) s5 = s6+r(i)*alpha(i,j)*a(i,j)**2*a(i-1,j)*br(i,j-B)*dr-8. #D0*alpha(i,j)*a(i,j)*c(i,j)**2*dr*dz*bz(i-1,j)-6.D0*r(i)*a #lpha(i,j)*a(i,j)*c(i,j)*dz*bz(i-1,j)*c(i+1,j)+6.D0*r(i)*alph #a(i,j)*a(i,j)*c(i,j)*dz*bz(i-1,j)*c(i-1,j)-3.D0*r(i)*alpha( #i,j)*c(i,j)**2*a(i-1,j)*dz*bz(i-1,j)-3.D0*r(i)*alpha(i,j)* #c(i,j)**2*a(i+1,j)*dz*bz(i+1,j) s6 = s5-r(i)*alpha(i,j)*a(i,j)**2*a(i+1,j)*br(i,j-B)*dr-16 #.D0*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i,j)+6.D0*r(i) #*alpha(i,j)*a(i,j)*c(i,j)*dz*bz(i+1,j)*c(i+1,j)-6.D0*r(i)*al #pha(i,j)*a(i,j)*c(i,j)*dz*bz(i+1,j)*c(i-1,j)+8.D0*r(i)*alpha #(i,j)*a(i,j)*c(i,j)**2*dz*bz(i+1,j) s4 = s6+8.D0*alpha(i,j)*a(i,j)*c(i,j)**2*dr*dz*bz(i+1,j)-8. #D0*alpha(i,j)*a(i,j)**2*br(i,j)*dr**2*a(i,j+W)+8.D0*alpha( #i,j)*a(i,j)**2*br(i,j)*dr**2*a(i,j-B)-2.D0*r(i)*a(i,j)*c( #i,j)**2*dz*alpha(i-1,j)*bz(i-1,j)+2.D0*r(i)*a(i,j)*c(i,j)* #*2.D0*dz*alpha(i-1,j)*bz(i+1,j)+2.D0*r(i)*a(i,j)*c(i,j)**2 #*dz*alpha(i+1,j)*bz(i-1,j)-2.D0*r(i)*a(i,j)*c(i,j)**2*dz*a #lpha(i+1,j)*bz(i+1,j) s5 = 1.D0/dr**2 s3 = s4*s5 s1 = s2*s3 s2 = 1.D0/dz/r(i)/alpha(i,j)**2/a(i,j)**4 Jac(2,4) = s1*s2 Jac(3,4) = -(-a(i,j)*alpha(i,j+W)+a(i,j)*alpha(i,j-B)-alpha(i,j) #*a(i,j-B)+alpha(i,j)*a(i,j+W))/a(i,j)/dz/r(i)/alpha(i,j)**2. #D0 Jac(4,4) = -1.D0/alpha(i,j)/a(i,j)**2*c(i,j)**2/dr**2 C---------------------------------------------------------------------- C Calculation of the residual. (Values/definitions of s1, ..., s6 C and elements of array res are stored in file residuals.) C---------------------------------------------------------------------- s2 = 1.D0/4.D0 s6 = c(i,j)**3*r(i)**2*dz**2*a(i-1,j)*alpha(i+1,j)+32.D0 #*a(i,j)**3*r(i)**2*alpha(i,j)*cpi*c(i,j)*phi_z(i,j)**2 #*dz**2*dr**2+32.D0*c(i,j)**3*r(i)**2*alpha(i,j)*cpi*a(i #,j)*phi_r(i,j)**2*dz**2*dr**2-2.D0*c(i,j)**2*r(i)**2*al #pha(i,j)*dz**2*a(i+1,j)*c(i+1,j)+8.D0*a(i,j)*r(i)*c(i,j)**2 #*alpha(i,j)*dr*dz**2*c(i+1,j)-8.D0*a(i,j)*r(i)*c(i,j)**2 #*alpha(i,j)*dr*dz**2*c(i-1,j)-8.D0*c(i,j)**3*r(i)**2*a( #i,j)*dz**2*alpha(i,j)-4.D0*a(i,j)*c(i,j)**3*r(i)*dr*dz**2*a #lpha(i-1,j)+c(i,j)**3*r(i)**2*dz**2*a(i+1,j)*alpha(i-1, #j)-c(i,j)**3*r(i)**2*dz**2*a(i-1,j)*alpha(i-1,j) s5 = s6+4.D0*c(i,j)**3*r(i)**2*a(i,j)*dz**2*alpha(i+1,j) #+8.D0*a(i,j)*c(i,j)**3*alpha(i,j)*dz**2*dr**2-c(i,j)**3* #r(i)**2*dz**2*a(i+1,j)*alpha(i+1,j)+a(i,j)*r(i)**2*c(i #,j)**2*dz**2*alpha(i+1,j)*c(i+1,j)+4.D0*a(i,j)*c(i,j)**3 #*r(i)*dr*dz**2*alpha(i+1,j)-a(i,j)*r(i)**2*c(i,j)**2*dz #**2*alpha(i+1,j)*c(i-1,j)-2.D0*c(i,j)**2*r(i)**2*alpha( #i,j)*dz**2*a(i-1,j)*c(i-1,j)+2.D0*r(i)**2*a(i,j)**2*alph #a(i,j)*dr**2*a(i,j+W)*c(i,j-B)+2.D0*r(i)**2*a(i,j)**2*a #lpha(i,j)*dr**2*a(i,j-B)*c(i,j+W)+8.D0*c(i,j)*r(i)**2*a(i #,j)**2*alpha(i,j)*a(i,j-B)*dr**2-2.D0*r(i)**2*a(i,j)**2 #*alpha(i,j)*dr**2*a(i,j+W)*c(i,j+W) s6 = s5-a(i,j)*r(i)**2*c(i,j)**2*dz**2*alpha(i-1,j)*c(i #+1,j)-24.D0*a(i,j)**3*r(i)**2*c(i,j)*dr**2*alpha(i,j)+a #(i,j)*r(i)**2*c(i,j)**2*dz**2*alpha(i-1,j)*c(i-1,j)+a(i #,j)**3*r(i)**2*dr**2*alpha(i,j+W)*c(i,j-B)+c(i,j)*r(i)* #*2.D0*a(i,j)**2*dr**2*a(i,j-B)*alpha(i,j-B)-c(i,j)*r(i)**2 #*a(i,j)**2*dr**2*a(i,j+W)*alpha(i,j-B)+a(i,j)**3*r(i #)**2*dr**2*alpha(i,j-B)*c(i,j+W)-a(i,j)**3*r(i)**2*dr #**2*alpha(i,j-B)*c(i,j-B)+4.D0*a(i,j)**3*r(i)**2*c(i,j) #*dr**2*alpha(i,j-B)+2.D0*c(i,j)**2*r(i)**2*alpha(i,j)*dz #**2*a(i+1,j)*c(i-1,j) s4 = s6+8.D0*c(i,j)*r(i)**2*a(i,j)**2*alpha(i,j)*a(i,j+ #W)*dr**2-8.D0*a(i,j)**3*c(i,j)**3*alpha(i,j)*dz**2*dr**2- #8.D0*c(i,j)**3*r(i)*alpha(i,j)*dz**2*a(i+1,j)*dr-2.D0*r(i) #**2*a(i,j)**2*alpha(i,j)*dr**2*a(i,j-B)*c(i,j-B)+8.D0*c #(i,j)**3*r(i)*alpha(i,j)*dz**2*a(i-1,j)*dr+2.D0*c(i,j)**2 #*r(i)**2*alpha(i,j)*dz**2*a(i-1,j)*c(i+1,j)+4.D0*c(i,j)** #3*r(i)**2*a(i,j)*dz**2*alpha(i-1,j)+c(i,j)*r(i)**2*a #(i,j)**2*dr**2*a(i,j+W)*alpha(i,j+W)+4.D0*a(i,j)**3*r(i #)**2*c(i,j)*dr**2*alpha(i,j+W)-c(i,j)*r(i)**2*a(i,j)**2 #*dr**2*a(i,j-B)*alpha(i,j+W)-a(i,j)**3*r(i)**2*dr**2* #alpha(i,j+W)*c(i,j+W) s5 = 1.D0/dr**2 s3 = s4*s5 s1 = s2*s3 s2 = 1.D0/a(i,j)**3/r(i)**2/c(i,j)**3/dz**2 res(1) = s1*s2 s4 = r(i)**2*c(i,j)**5*a(i,j)**2*dz**2*bz(i+1,j)*bz( #i-1,j)/4.D0+r(i)**2*c(i,j)**3*a(i,j)**4*dz*br(i,j-B)* #dr*bz(i+1,j)/4.D0-r(i)**2*c(i,j)**3*a(i,j)**4*dz*br( #i,j+W)*dr*bz(i+1,j)/4.D0-r(i)**2*c(i,j)*a(i,j)**6*dr**2* #br(i,j-B)**2/8.D0-2.D0*a(i,j)**3*c(i,j)*r(i)**2*alpha #(i,j)**2*a(i,j-B)*dr**2+a(i,j)**3*r(i)**2*alpha(i,j) #**2*dr**2*a(i,j-B)*c(i,j-B)/2.D0-c(i,j)**2*r(i)**2*al #pha(i,j)**2*a(i,j)*dz**2*a(i-1,j)*c(i+1,j)/2.D0+c(i,j)**2 #*r(i)**2*alpha(i,j)**2*a(i,j)*dz**2*a(i-1,j)*c(i-1,j) #/2.D0 s3 = s4-2.D0*c(i,j)**3*r(i)*alpha(i,j)**2*a(i,j)*dz**2*a #(i-1,j)*dr+r(i)**2*c(i,j)*a(i,j)**6*dr**2*br(i,j+W)*br #(i,j-B)/4.D0+r(i)**2*c(i,j)**3*a(i,j)**4*dz*br(i,j+W) #*dr*bz(i-1,j)/4.D0-r(i)**2*c(i,j)**5*a(i,j)**2*dz**2* #bz(i+1,j)**2/8.D0-r(i)**2*c(i,j)**3*a(i,j)**4*dz*b #r(i,j-B)*dr*bz(i-1,j)/4.D0-8.D0*cpi*a(i,j)**4*r(i)**2*c #(i,j)**3*alpha(i,j)**2*pi(i,j)**2*dz**2*dr**2-r(i)**2 #*c(i,j)**5*a(i,j)**2*dz**2*bz(i-1,j)**2/8.D0-8.D0*cpi #*a(i,j)**2*r(i)**2*c(i,j)**3*alpha(i,j)**2*phi_r(i #,j)**2*dz**2*dr**2 s4 = s3-8.D0*cpi*a(i,j)**4*r(i)**2*c(i,j)*alpha(i,j)**2 #*phi_z(i,j)**2*dz**2*dr**2-r(i)**2*c(i,j)*a(i,j)**6* #dr**2*br(i,j+W)**2/8.D0-4.D0*r(i)*krr(i,j)*c(i,j)**3*br #(i,j)*a(i,j)**2*alpha(i,j)*dz**2*dr**2-6.D0*br(i,j)**2*a #(i,j)**4*c(i,j)**3*dz**2*dr**2-2.D0*r(i)**2*krr(i,j)** #2.D0*c(i,j)**3*alpha(i,j)**2*dz**2*dr**2+2.D0*a(i,j)**4 #*c(i,j)**3*alpha(i,j)**2*dz**2*dr**2-2.D0*a(i,j)**2*c #(i,j)**3*alpha(i,j)**2*dz**2*dr**2-c(i,j)**2*r(i)**2 #*alpha(i,j)**2*a(i,j)*dz**2*a(i+1,j)*c(i-1,j)/2.D0 s2 = s4-2.D0*a(i,j)**2*c(i,j)**2*r(i)*dr*alpha(i,j)**2 #*dz**2*c(i+1,j)+2.D0*c(i,j)**3*r(i)*alpha(i,j)**2*a(i, #j)*dz**2*a(i+1,j)*dr+c(i,j)**2*r(i)**2*alpha(i,j)**2* #a(i,j)*dz**2*a(i+1,j)*c(i+1,j)/2.D0+2.D0*a(i,j)**2*c(i,j)* #*2.D0*r(i)*dr*alpha(i,j)**2*dz**2*c(i-1,j)-2.D0*a(i,j)**3 #*c(i,j)*r(i)**2*alpha(i,j)**2*a(i,j+W)*dr**2-a(i,j)**3 #*r(i)**2*alpha(i,j)**2*dr**2*a(i,j+W)*c(i,j-B)/2.D0-a #(i,j)**3*r(i)**2*alpha(i,j)**2*dr**2*a(i,j-B)*c(i,j+ #W)/2.D0+a(i,j)**3*r(i)**2*alpha(i,j)**2*dr**2*a(i,j+W #)*c(i,j+W)/2.D0+4.D0*a(i,j)**4*c(i,j)*r(i)**2*alpha(i,j #)**2*dr**2 s3 = 1.D0/dr**2/dz**2 s1 = s2*s3 s2 = 1.D0/a(i,j)**4/r(i)**2/c(i,j)**3/alpha(i,j)**2.D #0 res(2) = s1*s2 s4 = r(i)**2*alpha(i,j)*a(i,j)**5*dr*br(i,j-B)*c(i,j+W) #/8.D0-r(i)**2*c(i,j)**3*a(i,j)**3*alpha(i,j)*dz*bz(i #+1,j-B)/8.D0+c(i,j)**2*alpha(i,j)*r(i)*a(i,j)**3*br(i,j #)*dz**2*c(i+1,j)-r(i)**2*c(i,j)*a(i,j)**5*alpha(i,j)*dr #*br(i,j)-c(i,j)**2*alpha(i,j)*r(i)*a(i,j)**3*br(i,j)*d #z**2*c(i-1,j)-2.D0*c(i,j)**3*alpha(i,j)*a(i,j)**3*br(i, #j)*dr*dz**2-3.D0/8.D0*r(i)**2*c(i,j)*alpha(i,j)*a(i,j)**4 #*dr*a(i,j+W)*br(i,j-B)-r(i)**2*alpha(i,j)*a(i,j)**5*dr #*br(i,j-B)*c(i,j-B)/8.D0+r(i)**2*alpha(i,j)*a(i,j)**5*d #r*br(i,j+W)*c(i,j-B)/8.D0-r(i)**2*c(i,j)*a(i,j)**5*dr*a #lpha(i,j+W)*br(i,j+W)/8.D0 s3 = s4+r(i)**2*c(i,j)*a(i,j)**5*alpha(i,j)*dr*br(i,j+W #)/2.D0+3.D0/8.D0*r(i)**2*c(i,j)*alpha(i,j)*a(i,j)**4*dr* #a(i,j+W)*br(i,j+W)-r(i)**2*alpha(i,j)*a(i,j)**5*dr*br( #i,j+W)*c(i,j+W)/8.D0-3.D0/8.D0*r(i)**2*c(i,j)*alpha(i,j)*a #(i,j)**4*dr*a(i,j-B)*br(i,j+W)+r(i)**2*c(i,j)*a(i,j)**5 #*dr*alpha(i,j+W)*br(i,j-B)/8.D0-r(i)**2*c(i,j)*a(i,j)** #5.D0*dr*alpha(i,j-B)*br(i,j-B)/8.D0+r(i)**2*c(i,j)*a(i,j)* #*5.D0*alpha(i,j)*dr*br(i,j-B)/2.D0+3.D0/8.D0*r(i)**2*c(i,j) #*alpha(i,j)*a(i,j)**4*dr*a(i,j-B)*br(i,j-B)+c(i,j)**3*r #(i)**2*alpha(i,j)**2*krr(i,j)*a(i-1,j)*dz**2+r(i)**2 #*c(i,j)*a(i,j)**5*dr*alpha(i,j-B)*br(i,j+W)/8.D0 s4 = s3+r(i)**2*krr(i,j)*c(i,j)**2*alpha(i,j)**2*a(i #,j)*dz**2*c(i+1,j)-r(i)**2*krr(i,j)*c(i,j)**2*alpha(i,j #)**2*a(i,j)*dz**2*c(i-1,j)+2.D0*c(i,j)**3*krr(i,j)*alpha #(i,j)**2*r(i)*a(i,j)*dr*dz**2-c(i,j)**3*r(i)**2*alph #a(i,j)**2*krr(i,j)*a(i+1,j)*dz**2+8.D0*phi_r(i,j)*pi(i,j)*cp #i*r(i)**2*alpha(i,j)**2*a(i,j)**3*c(i,j)**3*dr*dz* #*2-r(i)**2*c(i,j)**3*alpha(i,j)*a(i,j)**2*a(i,j+W)*d #z*bz(i-1,j)/8.D0-r(i)**2*alpha(i,j)*dz*a(i,j)**3*c(i,j) #**2*bz(i+1,j)*c(i,j-B)/8.D0+r(i)**2*alpha(i,j)*dz*a(i,j #)**3*c(i,j)**2*bz(i+1,j)*c(i,j+W)/8.D0+r(i)**2*c(i,j #)**3*a(i,j)**3*alpha(i,j)*dz*bz(i-1,j-B)/8.D0 s5 = s4-r(i)**2*c(i,j)**3*a(i,j)**3*alpha(i,j-B)*dz*b #z(i-1,j)/8.D0-r(i)**2*alpha(i,j)*dz*a(i,j)**3*c(i,j)**2 #*bz(i-1,j)*c(i,j+W)/8.D0+r(i)**2*c(i,j)**3*alpha(i,j #)*a(i,j)**2*a(i,j-B)*dz*bz(i-1,j)/8.D0+r(i)**2*c(i,j)** #3*a(i,j)**3*alpha(i,j+W)*dz*bz(i-1,j)/8.D0+r(i)**2*al #pha(i,j)*dz*a(i,j)**3*c(i,j)**2*bz(i-1,j)*c(i,j-B)/8.D0 s2 = s5-r(i)**2*c(i,j)**3*a(i,j)**3*alpha(i,j)*dz*bz #(i-1,j+W)/8.D0+r(i)**2*c(i,j)**3*a(i,j)**3*alpha(i,j) #*dz*bz(i+1,j+W)/8.D0-r(i)**2*c(i,j)**3*alpha(i,j)*a(i,j #)**2*a(i,j-B)*dz*bz(i+1,j)/8.D0-r(i)**2*c(i,j)**3*a( #i,j)**3*alpha(i,j+W)*dz*bz(i+1,j)/8.D0+r(i)**2*c(i,j)**3 #*alpha(i,j)*a(i,j)**2*a(i,j+W)*dz*bz(i+1,j)/8.D0+r(i)** #2*c(i,j)**3*a(i,j)**3*alpha(i,j-B)*dz*bz(i+1,j)/8.D0 s3 = 1.D0/r(i)**2/alpha(i,j)**2 s1 = s2*s3 s2 = 1.D0/a(i,j)**3/c(i,j)**3/dr/dz**2 res(3) = s1*s2 s4 = -r(i)*a(i,j)*c(i,j)**3*dz*alpha(i-1,j)*bz(i-1,j)/8.D0 #+r(i)*alpha(i,j)*c(i,j)**3*a(i+1,j)*dz*bz(i-1,j)/8.D0+8.D0 #*phi_z(i,j)*pi(i,j)*cpi*c(i,j)*r(i)*alpha(i,j)**2*a(i,j)**3 #*dr**2*dz-r(i)*alpha(i,j)*c(i,j)**3*a(i-1,j)*dz*bz(i-1, #j)/8.D0+r(i)*a(i,j)*c(i,j)**3*dz*alpha(i+1,j)*bz(i-1,j)/8. #D0+r(i)*alpha(i,j)*a(i,j)*c(i,j)**3*dz*bz(i-1,j)/2.D0-alph #a(i,j)*a(i,j)*c(i,j)**3*dr*dz*bz(i-1,j)/2.D0+3.D0/8.D0*r(i #)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i-1,j)*c(i-1,j)-3.D0/8 #.D0*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i-1,j)*c(i+1,j #)-a(i,j)**3*c(i,j)*br(i,j)*dr**2*alpha(i,j-B) s3 = s4-r(i)*dr*alpha(i,j)*a(i,j)**3*br(i,j+W)*c(i-1,j)/8. #D0+r(i)*dr*alpha(i,j)*a(i,j)**3*br(i,j+W)*c(i+1,j)/8.D0+r #(i)*a(i,j)*c(i,j)**3*dz*alpha(i-1,j)*bz(i+1,j)/8.D0+r(i)*c #(i,j)*alpha(i,j)*a(i,j)**2*a(i+1,j)*br(i,j+W)*dr/8.D0-r(i #)*a(i,j)**3*c(i,j)*alpha(i+1,j)*br(i,j+W)*dr/8.D0+r(i)*a( #i,j)**3*c(i,j)*alpha(i-1,j)*br(i,j+W)*dr/8.D0-r(i)*c(i,j)* #alpha(i,j)*a(i,j)**2*a(i-1,j)*br(i,j+W)*dr/8.D0+alpha(i,j) #*a(i,j)*c(i,j)**3*dr*dz*bz(i+1,j)/2.D0-alpha(i,j)*a(i,j)** #3*c(i,j)*dr**2*br(i,j+W)/2.D0+r(i)*alpha(i,j)*a(i,j)**3 #*c(i,j)*br(i-1,j-B)*dr/8.D0 s4 = s3+r(i)*dr*alpha(i,j)*a(i,j)**3*br(i,j-B)*c(i-1,j)/8. #D0-r(i)*dr*alpha(i,j)*a(i,j)**3*br(i,j-B)*c(i+1,j)/8.D0+r #(i)*a(i,j)**3*c(i,j)*alpha(i+1,j)*br(i,j-B)*dr/8.D0-r(i)*a #lpha(i,j)*a(i,j)**3*c(i,j)*br(i-1,j+W)*dr/8.D0-r(i)*alpha #(i,j)*a(i,j)**3*c(i,j)*br(i+1,j-B)*dr/8.D0-c(i,j)*alpha(i, #j)*a(i,j)**2*br(i,j)*dr**2*a(i,j+W)+c(i,j)*alpha(i,j)*a(i #,j)**2*br(i,j)*dr**2*a(i,j-B)+r(i)*alpha(i,j)*a(i,j)**3 #*c(i,j)*br(i+1,j+W)*dr/8.D0+a(i,j)**3*c(i,j)*br(i,j)*dr**2 #*alpha(i,j+W) s5 = s4+alpha(i,j)*a(i,j)**3*c(i,j)*dr**2*br(i,j-B)/2.D0-r #(i)*a(i,j)**3*c(i,j)*alpha(i-1,j)*br(i,j-B)*dr/8.D0-r(i)*a #lpha(i,j)*a(i,j)*c(i,j)**3*dz*bz(i,j)+r(i)*c(i,j)*alpha( #i,j)*a(i,j)**2*a(i-1,j)*br(i,j-B)*dr/8.D0-r(i)*c(i,j)*alph #a(i,j)*a(i,j)**2*a(i+1,j)*br(i,j-B)*dr/8.D0 s2 = s5+r(i)*alpha(i,j)*c(i,j)**3*a(i-1,j)*dz*bz(i+1,j)/8. #D0-r(i)*alpha(i,j)*c(i,j)**3*a(i+1,j)*dz*bz(i+1,j)/8.D0-3. #D0/8.D0*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i+1,j)*c(i #-1,j)-r(i)*a(i,j)*c(i,j)**3*dz*alpha(i+1,j)*bz(i+1,j)/8.D0 #+r(i)*alpha(i,j)*a(i,j)*c(i,j)**3*dz*bz(i+1,j)/2.D0+3.D0/8 #.D0*r(i)*alpha(i,j)*a(i,j)*c(i,j)**2*dz*bz(i+1,j)*c(i+1,j #) s3 = 1.D0/dr**2/dz s1 = s2*s3 s2 = 1.D0/c(i,j)/r(i)/alpha(i,j)**2/a(i,j)**3 res(4) = s1*s2 C---------------------------------------------------------------------- C Tracing statement. C---------------------------------------------------------------------- if ( ltrace1 ) then write(0,*) 'res(',i,',',j,') =',res endif res(1) = res(1) - mg_rhs1(i,j) res(2) = res(2) - mg_rhs2(i,j) res(3) = res(3) - mg_rhs3(i,j) res(4) = res(4) - mg_rhs4(i,j) C---------------------------------------------------------------------- C The linear solver solves Jac * du = res. C---------------------------------------------------------------------- call linear_solver(Jac,res,du,4,info) if (info .ne. 0) then write(0,*) 'i=',i,', j=',j stop endif C---------------------------------------------------------------------- C Update the values of the function at the point. C---------------------------------------------------------------------- alpha(i,j) = alpha(i,j) - du(1) a(i,j) = a(i,j) - du(2) br(i,j) = br(i,j) - du(3) bz(i,j) = bz(i,j) - du(4) do k = 1, 4 res_relax = res_relax + res(k)**2 / 4.0d0**2 enddo if ( ltrace0 ) then write(0,*) cdnm,'=', res_relax, & ', (i,j)=(',i,',',j,')' write(0,*) 'du=',du write(0,*) 'Sources=',phi(i,j),c(i,j),krr(i,j),pi(i,j) endif end do if (order .eq. 1) isw = 3 -isw end do if (order .eq. 1) jsw = 3 - jsw end do C---------------------------------------------------------------------- C Neumman Boundary Conditions: C C alpha @ r=0 C C 1/r fall off i.e. f'+f/r=0 C C a @ r=rmax C br @ r=rmax C bz @ r=rmax C C Special for boundary condition for alpha: C C alpha = 1/a @ r=rmax C---------------------------------------------------------------------- C---------------------------------------------------------------------- C Set the boundary conditions at r=rmax for alpha. C---------------------------------------------------------------------- C do j = 1 , Nz C alpha(Nr,j) = 1.0d0 / a(Nr,j) C end do do j = 1 , Nz C a(Nr, j) = (4.0d0 * a(Nr-1,j) - a(Nr-2,j)) / 3.0d0 a(Nr,j) = (2.0d0 * a(Nr-1,j) * r(Nr) / dr & - 0.5d0 * a(Nr-2,j) * r(Nr) / dr + 1.0d0 ) / & (1.5d0 * r(Nr) / dr + 1.0d0 ) br(Nr,j) = (2.0d0 * br(Nr-1,j) / dr - 0.5d0 * br(Nr-2,j)/dr)/ & (1.5d0 / dr + 1.0d0 / r(Nr)) bz(Nr,j) = (2.0d0 * bz(Nr-1,j) / dr - 0.5d0 * bz(Nr-2,j)/dr)/ & (1.5d0 / dr + 1.0d0 / r(Nr)) alpha(1, j) = (4.0d0 * alpha(2,j) - alpha(3,j)) / 3.0d0 alpha(Nr,j) = (2.0d0 * alpha(Nr-1,j) * r(Nr) / dr & - 0.5d0 * alpha(Nr-2,j) * r(Nr) / dr + 1.0d0 ) / & (1.5d0 * r(Nr) / dr + 1.0d0 ) enddo do i = 2, Nr do j = 1 , Nz krr(i,j) = -2.0d0 * a(i,j)**2 * br(i,j) & / (alpha(i,j) * r(i)) end do end do res_relax = sqrt(res_relax / (Nr * Nz)) return C---------------------------------------------------------------------- C End function. C---------------------------------------------------------------------- end