        subroutine  calc_resid(alpha, br, bz, a, c, krr, 
     &                         phi, pi, phi_r, phi_z, r, z,
     &                         Nr, Nz, res1, res2, res3, res4)

        implicit      none 

        logical       ltrace
        parameter   ( ltrace = .true. )  
          
        character*10  cdnm 
        parameter   ( cdnm = 'calc_resid' ) 

        integer       Nr,    Nz
         
        real*8        r(nr),    z(Nz)

        real*8        alpha(Nr, Nz)
        real*8        br(Nr, Nz)
        real*8        bz(Nr, Nz)
        real*8        a(Nr, Nz)
        real*8        c(Nr, Nz)
        real*8        krr(Nr, Nz)

        real*8        phi(Nr, Nz)
        real*8        pi(Nr, Nz)
        real*8        phi_r(Nr, Nz)
        real*8        phi_z(Nr, Nz)

        real*8        s1
        real*8        s2
        real*8        s3
        real*8        s4
        real*8        s5
        real*8        s6

        real*8        res1(Nr, Nz)         
        real*8        res2(Nr, Nz)         
        real*8        res3(Nr, Nz)         
        real*8        res4(Nr, Nz)         

        integer       i,     j       

        integer       W,     B

        real*8        dr,    dz
 
        real*8        cpi
        parameter   ( cpi = 3.141 5926 5358 9793 d0 )

        real*8        norm_res1
        real*8        norm_res2
        real*8        norm_res3
        real*8        norm_res4

c===================================================================
c===================================================================
      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
      endif
      norm_res1 = 0.0d0 
      norm_res2 = 0.0d0
      norm_res3 = 0.0d0
      norm_res4 = 0.0d0
      do j = 1, Nz   
         do i = 2, Nr-1
            W = 1
            B = 1
            if ( j .eq. 1)   B = j - Nz             
            if ( j .eq. Nz)  W = -j + 1
      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.D
     #0*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.
     #D0*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.
     #D0*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
     #.D0*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.D
     #0*r(i)**2*alpha(i,j)*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*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
     #.D0*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
      res1(i,j) = 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.
     #D0*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.D
     #0*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.
     #D0*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.D
     #0*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.
     #D0*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.D
     #0*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.D
     #0*c(i,j)*r(i)**2*alpha(i,j)**2*a(i,j+W)*dr**2-a(i,j)**3
     #.D0*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
      res2(i,j) = 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.D
     #0*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
     #.D0*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
     #.D0*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.D0*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
     #.D0*alpha(i,j)*a(i,j)**2*a(i,j+W)*dz*bz(i+1,j)/8.D0+r(i)**
     #2.D0*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
      res3(i,j) = 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
     #.D0*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.D0*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
      res4(i,j) = s1*s2
             norm_res1 = norm_res1 + res1(i,j)**2 
             norm_res2 = norm_res2 + res2(i,j)**2 
             norm_res3 = norm_res3 + res3(i,j)**2 
             norm_res4 = norm_res4 + res4(i,j)**2 
         enddo
      enddo 
      norm_res1 = sqrt(norm_res1/(Nr*Nz)) 
      norm_res2 = sqrt(norm_res2/(Nr*Nz)) 
      norm_res3 = sqrt(norm_res3/(Nr*Nz)) 
      norm_res4 = sqrt(norm_res4/(Nr*Nz)) 
      if ( ltrace ) then
         write(0,*) cdnm,': norm_res1=',norm_res1
         write(0,*) cdnm,': norm_res2=',norm_res2
         write(0,*) cdnm,': norm_res3=',norm_res3
         write(0,*) cdnm,': norm_res4=',norm_res4
      endif
      return
      end
