[comp.sys.apollo] FORTRAN COMPILER BUG

rogden@uceng.UC.EDU (rob ogden) (03/20/90)

Hello  comp.sys.apollo,

FORTRAN compiler bug  DN10000:

Here is a FORTRAN bug that was found on out dn10000 last week.
There are ~560 lines in this email because the problem subroutine is
attached. I heard this morning, 3/19/90, from Debbie Sahrbeck 
of Apollo.
The work around is:

compile with    -save
with the  f77  ,unix command , this needs to be prefixed with  -W0
example:  f77 -W0,-save -c chrstf.f

An Apollo Problem Report number has been assigned. APR#  DE0CE
   -      -       -

We have a DN10000 with FORTRAN v10.7.p
Debbie Sahrbeck indicates that the problem does NOT show up on
the 68000s.

That which follows is the information sent to Debbie Sahrbeck.

Rob Ogden, rogden@uceng.uc.edu
           uccba!uceng!rogden
--------------------------------------------------------------

Hello Debbie Sahrbeck, sahrbeck_d@apollo.com

I am Rob Ogden, rogden@uceng.uc.edu.
The following information pertains to reference #A2015160
a FORTRAN compiler error.

Attached is the subroutine chrstf.
The compiler errors are shown.

The following is a table of parameters .

This works          ->  parameter (ngrd=1,imx=13,jmx=63,kmx=9)
This does not work  ->  parameter (ngrd=1,imx=13,jmx=64,kmx=9)
This does not work  ->  parameter (ngrd=1,imx=13,jmx=63,kmx=10)

nasa.hm.ztst3d[80] f77 -c chrstf.f
chrstf.f:
?(ftn) LIR_$BUILD_OPERAND - signed_immed16 out of range: 4294907296, ADD
Compiler backend failure while processing routine "chrstf" 
?(ftn) All routines processed by backend - aborting due to previous errors
(00594)       end
**** Error #319 on Line 594: System fault trapped, fault message follows:
status 7F000000
Apollo-specific fault
Compiler error in file chrstf.f: see a system manager
nasa.hm.ztst3d[81] 
nasa.hm.ztst3d[81] 
nasa.hm.ztst3d[81] f77 -O -c chrstf.f
chrstf.f:
**** Informational #929 on Line 35: Assignment eliminated, value never used: ybody
**** Informational #929 on Line 29: Assignment eliminated, value never used: dfdx
**** Informational #929 on Line 21: Assignment eliminated, value never used: dfdx
**** Informational #929 on Line 20: Assignment eliminated, value never used: fofx
**** Informational #929 on Line 19: Assignment eliminated, value never used: ybody
**** Informational #929 on Line 17: Assignment eliminated, value never used: yny
?(ftn) LIR_$BUILD_OPERAND - signed_immed16 out of range: 4294907296, ADD
Compiler backend failure while processing routine "chrstf" 
?(ftn) All routines processed by backend - aborting due to previous errors
(00594)       end
**** Error #319 System fault trapped, fault message follows:
status 7F000000
Apollo-specific fault
Compiler error in file chrstf.f: see a system manager
nasa.hm.ztst3d[82] 
nasa.hm.ztst3d[82] 
nasa.hm.ztst3d[82] /com/sh
$ 
$ ftn chrstf.ftn
**** Informational #929 on Line 35: Assignment eliminated, value never used: ybody
**** Informational #929 on Line 29: Assignment eliminated, value never used: dfdx
**** Informational #929 on Line 21: Assignment eliminated, value never used: dfdx
**** Informational #929 on Line 20: Assignment eliminated, value never used: fofx
**** Informational #929 on Line 19: Assignment eliminated, value never used: ybody
**** Informational #929 on Line 17: Assignment eliminated, value never used: yny
?(ftn) LIR_$BUILD_OPERAND - signed_immed16 out of range: 4294907296, ADD
Compiler backend failure while processing routine "chrstf" 
?(ftn) All routines processed by backend - aborting due to previous errors
(00594)       end
**** Error #319 System fault trapped, fault message follows:
status 7F000000
?(sh) "/com/ftn" - status 7F000000
$ 

THE SUBROUTINE STARTS HERE.
---------------------------
please note the parameter statement, the third line in the subroutine.
The following is a table of parameter statements that have been tried.

This works          ->  parameter (ngrd=1,imx=13,jmx=63,kmx=9)
This does not work  ->  parameter (ngrd=1,imx=13,jmx=64,kmx=9)
This does not work  ->  parameter (ngrd=1,imx=13,jmx=63,kmx=10)
-----------------------------------

      subroutine chrstf (ip,i)
c      implicit real*8(a-h,o-z)
      parameter (ngrd=1,imx=13,jmx=64,kmx=9)
      common /grsz/ nx,ny,nz,nw,ns
      common /grid/ x(imx),y(0:jmx+1),z(0:kmx+1)
     1             ,hx(imx),hy(0:jmx+1),hz(0:kmx+1)
      common /coords/ xp(0:jmx+1,0:kmx+1,3),yp(0:jmx+1,0:kmx+1,3)
     1               ,zp(0:jmx+1,0:kmx+1,3)
      common /metten/ sqrg(jmx,kmx,2),gv11(jmx,kmx),gv12(jmx,kmx)
     1   ,gv13(jmx,kmx),gv22(jmx,kmx),gv23(jmx,kmx)
     2   ,gv33(jmx,kmx),gn11(jmx,kmx,2),gn12(jmx,kmx,2)
     3   ,gn13(jmx,kmx,2),gn22(jmx,kmx,2),gn23(jmx,kmx,2)
     4   ,gn33(jmx,kmx,2),chf(jmx,kmx,3,3,3)
      dimension dgdx(jmx,kmx,3,3,3)
c     
      dx=hx(ip)+hx(ip+1)
      do 90 k=2,nz-1
      dz=hz(k)+hz(k+1)
      do 90 j=2,ny-1
         dy=hy(j)+hy(j+1)
         dxpdx=(xp(j,k,i+1)-xp(j,k,i))/hx(ip+1)
         dypdx=(yp(j,k,i+1)-yp(j,k,i))/hx(ip+1)
         dzpdx=(zp(j,k,i+1)-zp(j,k,i))/hx(ip+1)
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k,i+1)-xp(j-1,k,i)
     1   -xp(j-1,k,i+1))/dy
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k,i+1)-yp(j-1,k,i)
     1   -yp(j-1,k,i+1))/dy
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k,i+1)-zp(j-1,k,i)
     1   -zp(j-1,k,i+1))/dy
         dxpdz=0.5e0*(xp(j,k+1,i)+xp(j,k+1,i+1)-xp(j,k-1,i)
     1   -xp(j,k-1,i+1))/dz
         dypdz=0.5e0*(yp(j,k+1,i)+yp(j,k+1,i+1)-yp(j,k-1,i)
     1   -yp(j,k-1,i+1))/dz
         dzpdz=0.5e0*(zp(j,k+1,i)+zp(j,k+1,i+1)-zp(j,k-1,i)
     1   -zp(j,k-1,i+1))/dz
         g11ip=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12ip=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13ip=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22ip=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23ip=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33ip=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=(xp(j,k,i)-xp(j,k,i-1))/hx(ip)
         dypdx=(yp(j,k,i)-yp(j,k,i-1))/hx(ip)
         dzpdx=(zp(j,k,i)-zp(j,k,i-1))/hx(ip)
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k,i-1)-xp(j-1,k,i)
     1   -xp(j-1,k,i-1))/dy
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k,i-1)-yp(j-1,k,i)
     1   -yp(j-1,k,i-1))/dy
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k,i-1)-zp(j-1,k,i)
     1   -zp(j-1,k,i-1))/dy
         dxpdz=0.5e0*(xp(j,k+1,i)+xp(j,k+1,i-1)-xp(j,k-1,i)
     1   -xp(j,k-1,i-1))/dz
         dypdz=0.5e0*(yp(j,k+1,i)+yp(j,k+1,i-1)-yp(j,k-1,i)
     1   -yp(j,k-1,i-1))/dz
         dzpdz=0.5e0*(zp(j,k+1,i)+zp(j,k+1,i-1)-zp(j,k-1,i)
     1   -zp(j,k-1,i-1))/dz
         g11im=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12im=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13im=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22im=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23im=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33im=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=0.5e0*(xp(j,k,i+1)+xp(j+1,k,i+1)-xp(j,k,i-1)
     1   -xp(j+1,k,i-1))/dx
         dypdx=0.5e0*(yp(j,k,i+1)+yp(j+1,k,i+1)-yp(j,k,i-1)
     1   -yp(j+1,k,i-1))/dx
         dzpdx=0.5e0*(zp(j,k,i+1)+zp(j+1,k,i+1)-zp(j,k,i-1)
     1   -zp(j+1,k,i-1))/dx
         dxpdy=(xp(j+1,k,i)-xp(j,k,i))/hy(j+1)
         dypdy=(yp(j+1,k,i)-yp(j,k,i))/hy(j+1)
         dzpdy=(zp(j+1,k,i)-zp(j,k,i))/hy(j+1)
         dxpdz=0.5e0*(xp(j,k+1,i)+xp(j+1,k+1,i)-xp(j,k-1,i)
     1   -xp(j+1,k-1,i))/dz
         dypdz=0.5e0*(yp(j,k+1,i)+yp(j+1,k+1,i)-yp(j,k-1,i)
     1   -yp(j+1,k-1,i))/dz
         dzpdz=0.5e0*(zp(j,k+1,i)+zp(j+1,k+1,i)-zp(j,k-1,i)
     1   -zp(j+1,k-1,i))/dz
         g11jp=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12jp=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13jp=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22jp=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23jp=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33jp=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=0.5e0*(xp(j,k,i+1)+xp(j-1,k,i+1)-xp(j,k,i-1)
     1   -xp(j-1,k,i-1))/dx
         dypdx=0.5e0*(yp(j,k,i+1)+yp(j-1,k,i+1)-yp(j,k,i-1)
     1   -yp(j-1,k,i-1))/dx
         dzpdx=0.5e0*(zp(j,k,i+1)+zp(j-1,k,i+1)-zp(j,k,i-1)
     1   -zp(j-1,k,i-1))/dx
         dxpdy=(xp(j,k,i)-xp(j-1,k,i))/hy(j)
         dypdy=(yp(j,k,i)-yp(j-1,k,i))/hy(j)
         dzpdy=(zp(j,k,i)-zp(j-1,k,i))/hy(j)
         dxpdz=0.5e0*(xp(j,k+1,i)+xp(j-1,k+1,i)-xp(j,k-1,i)
     1   -xp(j-1,k-1,i))/dz
         dypdz=0.5e0*(yp(j,k+1,i)+yp(j-1,k+1,i)-yp(j,k-1,i)
     1   -yp(j-1,k-1,i))/dz
         dzpdz=0.5e0*(zp(j,k+1,i)+zp(j-1,k+1,i)-zp(j,k-1,i)
     1   -zp(j-1,k-1,i))/dz
         g11jm=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12jm=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13jm=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22jm=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23jm=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33jm=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=0.5e0*(xp(j,k,i+1)+xp(j,k+1,i+1)-xp(j,k,i-1)
     1   -xp(j,k+1,i-1))/dx
         dypdx=0.5e0*(yp(j,k,i+1)+yp(j,k+1,i+1)-yp(j,k,i-1)
     1   -yp(j,k+1,i-1))/dx
         dzpdx=0.5e0*(zp(j,k,i+1)+zp(j,k+1,i+1)-zp(j,k,i-1)
     1   -zp(j,k+1,i-1))/dx
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k+1,i)-xp(j-1,k,i)
     1   -xp(j-1,k+1,i))/dy
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k+1,i)-yp(j-1,k,i)
     1   -yp(j-1,k+1,i))/dy
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k+1,i)-zp(j-1,k,i)
     1   -zp(j-1,k+1,i))/dy
         dxpdz=(xp(j,k+1,i)-xp(j,k,i))/hz(k+1)
         dypdz=(yp(j,k+1,i)-yp(j,k,i))/hz(k+1)
         dzpdz=(zp(j,k+1,i)-zp(j,k,i))/hz(k+1)
         g11kp=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12kp=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13kp=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22kp=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23kp=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33kp=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=0.5e0*(xp(j,k,i+1)+xp(j,k-1,i+1)-xp(j,k,i-1)
     1   -xp(j,k-1,i-1))/dx
         dypdx=0.5e0*(yp(j,k,i+1)+yp(j,k-1,i+1)-yp(j,k,i-1)
     1   -yp(j,k-1,i-1))/dx
         dzpdx=0.5e0*(zp(j,k,i+1)+zp(j,k-1,i+1)-zp(j,k,i-1)
     1   -zp(j,k-1,i-1))/dx
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k-1,i)-xp(j-1,k,i)
     1   -xp(j-1,k-1,i))/dy
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k-1,i)-yp(j-1,k,i)
     1   -yp(j-1,k-1,i))/dy
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k-1,i)-zp(j-1,k,i)
     1   -zp(j-1,k-1,i))/dy
         dxpdz=(xp(j,k,i)-xp(j,k-1,i))/hz(k)
         dypdz=(yp(j,k,i)-yp(j,k-1,i))/hz(k)
         dzpdz=(zp(j,k,i)-zp(j,k-1,i))/hz(k)
         g11km=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12km=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13km=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22km=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23km=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33km=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dgdx(j,k,1,1,1)=(g11ip-g11im)*2.0e0/dx
         dgdx(j,k,1,2,1)=(g12ip-g12im)*2.0e0/dx
         dgdx(j,k,1,3,1)=(g13ip-g13im)*2.0e0/dx
         dgdx(j,k,2,2,1)=(g22ip-g22im)*2.0e0/dx
         dgdx(j,k,2,3,1)=(g23ip-g23im)*2.0e0/dx
         dgdx(j,k,3,3,1)=(g33ip-g33im)*2.0e0/dx
         dgdx(j,k,2,1,1)=dgdx(j,k,1,2,1)
         dgdx(j,k,3,1,1)=dgdx(j,k,1,3,1)
         dgdx(j,k,3,2,1)=dgdx(j,k,2,3,1)
         dgdx(j,k,1,1,2)=(g11jp-g11jm)*2.0e0/dy
         dgdx(j,k,1,2,2)=(g12jp-g12jm)*2.0e0/dy
         dgdx(j,k,1,3,2)=(g13jp-g13jm)*2.0e0/dy
         dgdx(j,k,2,2,2)=(g22jp-g22jm)*2.0e0/dy
         dgdx(j,k,2,3,2)=(g23jp-g23jm)*2.0e0/dy
         dgdx(j,k,3,3,2)=(g33jp-g33jm)*2.0e0/dy
         dgdx(j,k,2,1,2)=dgdx(j,k,1,2,2)
         dgdx(j,k,3,1,2)=dgdx(j,k,1,3,2)
         dgdx(j,k,3,2,2)=dgdx(j,k,2,3,2)
         dgdx(j,k,1,1,3)=(g11kp-g11km)*2.0e0/dz
         dgdx(j,k,1,2,3)=(g12kp-g12km)*2.0e0/dz
         dgdx(j,k,1,3,3)=(g13kp-g13km)*2.0e0/dz
         dgdx(j,k,2,2,3)=(g22kp-g22km)*2.0e0/dz
         dgdx(j,k,2,3,3)=(g23kp-g23km)*2.0e0/dz
         dgdx(j,k,3,3,3)=(g33kp-g33km)*2.0e0/dz
         dgdx(j,k,2,1,3)=dgdx(j,k,1,2,3)
         dgdx(j,k,3,1,3)=dgdx(j,k,1,3,3)
         dgdx(j,k,3,2,3)=dgdx(j,k,2,3,3)
 90   continue
c     
      do 120 k=2,nz-1
      do 120 j=2,ny-1
         chf(j,k,1,1,1)=0.5e0*gn11(j,k,i)*
     1   (dgdx(j,k,1,1,1)+dgdx(j,k,1,1,1)-dgdx(j,k,1,1,1))
     2   +0.5e0*gn12(j,k,i)*
     3   (dgdx(j,k,1,2,1)+dgdx(j,k,2,1,1)-dgdx(j,k,1,1,2))
     4   +0.5e0*gn13(j,k,i)*
     5   (dgdx(j,k,1,3,1)+dgdx(j,k,3,1,1)-dgdx(j,k,1,1,3))
         chf(j,k,1,1,2)=0.5e0*gn11(j,k,i)*
     1   (dgdx(j,k,1,1,2)+dgdx(j,k,1,2,1)-dgdx(j,k,1,2,1))
     2   +0.5e0*gn12(j,k,i)*
     3   (dgdx(j,k,1,2,2)+dgdx(j,k,2,2,1)-dgdx(j,k,1,2,2))
     4   +0.5e0*gn13(j,k,i)*
     5   (dgdx(j,k,1,3,2)+dgdx(j,k,3,2,1)-dgdx(j,k,1,2,3))
         chf(j,k,1,1,3)=0.5e0*gn11(j,k,i)*
     1   (dgdx(j,k,1,1,3)+dgdx(j,k,1,3,1)-dgdx(j,k,1,3,1))
     2   +0.5e0*gn12(j,k,i)*
     3   (dgdx(j,k,1,2,3)+dgdx(j,k,2,3,1)-dgdx(j,k,1,3,2))
     4   +0.5e0*gn13(j,k,i)*
     5   (dgdx(j,k,1,3,3)+dgdx(j,k,3,3,1)-dgdx(j,k,1,3,3))
         chf(j,k,1,2,2)=0.5e0*gn11(j,k,i)*
     1   (dgdx(j,k,2,1,2)+dgdx(j,k,1,2,2)-dgdx(j,k,2,2,1))
     2   +0.5e0*gn12(j,k,i)*
     3   (dgdx(j,k,2,2,2)+dgdx(j,k,2,2,2)-dgdx(j,k,2,2,2))
     4   +0.5e0*gn13(j,k,i)*
     5   (dgdx(j,k,2,3,2)+dgdx(j,k,3,2,2)-dgdx(j,k,2,2,3))
         chf(j,k,1,2,3)=0.5e0*gn11(j,k,i)*
     1   (dgdx(j,k,2,1,3)+dgdx(j,k,1,3,2)-dgdx(j,k,2,3,1))
     2   +0.5e0*gn12(j,k,i)*
     3   (dgdx(j,k,2,2,3)+dgdx(j,k,2,3,2)-dgdx(j,k,2,3,2))
     4   +0.5e0*gn13(j,k,i)*
     5   (dgdx(j,k,2,3,3)+dgdx(j,k,3,3,2)-dgdx(j,k,2,3,3))
         chf(j,k,1,3,3)=0.5e0*gn11(j,k,i)*
     1   (dgdx(j,k,3,1,3)+dgdx(j,k,1,3,3)-dgdx(j,k,3,3,1))
     2   +0.5e0*gn12(j,k,i)*
     3   (dgdx(j,k,3,2,3)+dgdx(j,k,2,3,3)-dgdx(j,k,3,3,2))
     4   +0.5e0*gn13(j,k,i)*
     5   (dgdx(j,k,3,3,3)+dgdx(j,k,3,3,3)-dgdx(j,k,3,3,3))
         chf(j,k,1,2,1)=chf(j,k,1,1,2)
         chf(j,k,1,3,1)=chf(j,k,1,1,3)
         chf(j,k,1,3,2)=chf(j,k,1,2,3)
         chf(j,k,3,1,1)=0.5e0*gn13(j,k,i)*
     1   (dgdx(j,k,1,1,1)+dgdx(j,k,1,1,1)-dgdx(j,k,1,1,1))
     2   +0.5e0*gn23(j,k,i)*
     3   (dgdx(j,k,1,2,1)+dgdx(j,k,2,1,1)-dgdx(j,k,1,1,2))
     4   +0.5e0*gn33(j,k,i)*
     5   (dgdx(j,k,1,3,1)+dgdx(j,k,3,1,1)-dgdx(j,k,1,1,3))
         chf(j,k,3,1,2)=0.5e0*gn13(j,k,i)*
     1   (dgdx(j,k,1,1,2)+dgdx(j,k,1,2,1)-dgdx(j,k,1,2,1))
     2   +0.5e0*gn23(j,k,i)*
     3   (dgdx(j,k,1,2,2)+dgdx(j,k,2,2,1)-dgdx(j,k,1,2,2))
     4   +0.5e0*gn33(j,k,i)*
     5   (dgdx(j,k,1,3,2)+dgdx(j,k,3,2,1)-dgdx(j,k,1,2,3))
         chf(j,k,3,1,3)=0.5e0*gn13(j,k,i)*
     1   (dgdx(j,k,1,1,3)+dgdx(j,k,1,3,1)-dgdx(j,k,1,3,1))
     2   +0.5e0*gn23(j,k,i)*
     3   (dgdx(j,k,1,2,3)+dgdx(j,k,2,3,1)-dgdx(j,k,1,3,2))
     4   +0.5e0*gn33(j,k,i)*
     5   (dgdx(j,k,1,3,3)+dgdx(j,k,3,3,1)-dgdx(j,k,1,3,3))
         chf(j,k,3,2,2)=0.5e0*gn13(j,k,i)*
     1   (dgdx(j,k,2,1,2)+dgdx(j,k,1,2,2)-dgdx(j,k,2,2,1))
     2   +0.5e0*gn23(j,k,i)*
     3   (dgdx(j,k,2,2,2)+dgdx(j,k,2,2,2)-dgdx(j,k,2,2,2))
     4   +0.5e0*gn33(j,k,i)*
     5   (dgdx(j,k,2,3,2)+dgdx(j,k,3,2,2)-dgdx(j,k,2,2,3))
         chf(j,k,3,2,3)=0.5e0*gn13(j,k,i)*
     1   (dgdx(j,k,2,1,3)+dgdx(j,k,1,3,2)-dgdx(j,k,2,3,1))
     2   +0.5e0*gn23(j,k,i)*
     3   (dgdx(j,k,2,2,3)+dgdx(j,k,2,3,2)-dgdx(j,k,2,3,2))
     4   +0.5e0*gn33(j,k,i)*
     5   (dgdx(j,k,2,3,3)+dgdx(j,k,3,3,2)-dgdx(j,k,2,3,3))
         chf(j,k,3,3,3)=0.5e0*gn13(j,k,i)*
     1   (dgdx(j,k,3,1,3)+dgdx(j,k,1,3,3)-dgdx(j,k,3,3,1))
     2   +0.5e0*gn23(j,k,i)*
     3   (dgdx(j,k,3,2,3)+dgdx(j,k,2,3,3)-dgdx(j,k,3,3,2))
     4   +0.5e0*gn33(j,k,i)*
     5   (dgdx(j,k,3,3,3)+dgdx(j,k,3,3,3)-dgdx(j,k,3,3,3))
         chf(j,k,3,2,1)=chf(j,k,3,1,2)
         chf(j,k,3,3,1)=chf(j,k,3,1,3)
         chf(j,k,3,3,2)=chf(j,k,3,2,3)
 120  continue
c     
      dx=hx(ip)+hx(ip+1)
      do 150 k=2,nz-1
      dz=hz(k)+hz(k+1)
      do 150 j=1,ny-1
         dxpdx=0.5e0*(xp(j,k,i+1)+xp(j+1,k,i+1)-xp(j,k,i)
     1   -xp(j+1,k,i))/hx(ip+1)
         dypdx=0.5e0*(yp(j,k,i+1)+yp(j+1,k,i+1)-yp(j,k,i)
     1   -yp(j+1,k,i))/hx(ip+1)
         dzpdx=0.5e0*(zp(j,k,i+1)+zp(j+1,k,i+1)-zp(j,k,i)
     1   -zp(j+1,k,i))/hx(ip+1)
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k,i+1)-xp(j,k,i)
     1   -xp(j,k,i+1))/hy(j+1)
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k,i+1)-yp(j,k,i)
     1   -yp(j,k,i+1))/hy(j+1)
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k,i+1)-zp(j,k,i)
     1   -zp(j,k,i+1))/hy(j+1)
         dxpdz=0.25e0*(xp(j,k+1,i)+xp(j+1,k+1,i)+xp(j,k+1,i+1)
     1   +xp(j+1,k+1,i+1)-xp(j,k-1,i)-xp(j+1,k-1,i)-xp(j,k-1,i+1)
     2   -xp(j+1,k-1,i+1))/dz
         dypdz=0.25e0*(yp(j,k+1,i)+yp(j+1,k+1,i)+yp(j,k+1,i+1)
     1   +yp(j+1,k+1,i+1)-yp(j,k-1,i)-yp(j+1,k-1,i)-yp(j,k-1,i+1)
     2   -yp(j+1,k-1,i+1))/dz
         dzpdz=0.25e0*(zp(j,k+1,i)+zp(j+1,k+1,i)+zp(j,k+1,i+1)
     1   +zp(j+1,k+1,i+1)-zp(j,k-1,i)-zp(j+1,k-1,i)-zp(j,k-1,i+1)
     2   -zp(j+1,k-1,i+1))/dz
         g11ip=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12ip=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13ip=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22ip=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23ip=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33ip=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=0.5e0*(xp(j,k,i)+xp(j+1,k,i)-xp(j,k,i-1)
     1   -xp(j+1,k,i-1))/hx(ip)
         dypdx=0.5e0*(yp(j,k,i)+yp(j+1,k,i)-yp(j,k,i-1)
     1   -yp(j+1,k,i-1))/hx(ip)
         dzpdx=0.5e0*(zp(j,k,i)+zp(j+1,k,i)-zp(j,k,i-1)
     1   -zp(j+1,k,i-1))/hx(ip)
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k,i-1)-xp(j,k,i)
     1   -xp(j,k,i-1))/hy(j+1)
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k,i-1)-yp(j,k,i)
     1   -yp(j,k,i-1))/hy(j+1)
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k,i-1)-zp(j,k,i)
     1   -zp(j,k,i-1))/hy(j+1)
         dxpdz=0.25e0*(xp(j,k+1,i)+xp(j+1,k+1,i)+xp(j,k+1,i-1)
     1   +xp(j+1,k+1,i-1)-xp(j,k-1,i)-xp(j+1,k-1,i)-xp(j,k-1,i-1)
     2   -xp(j+1,k-1,i-1))/dz
         dypdz=0.25e0*(yp(j,k+1,i)+yp(j+1,k+1,i)+yp(j,k+1,i-1)
     1   +yp(j+1,k+1,i-1)-yp(j,k-1,i)-yp(j+1,k-1,i)-yp(j,k-1,i-1)
     2   -yp(j+1,k-1,i-1))/dz
         dzpdz=0.25e0*(zp(j,k+1,i)+zp(j+1,k+1,i)+zp(j,k+1,i-1)
     1   +zp(j+1,k+1,i-1)-zp(j,k-1,i)-zp(j+1,k-1,i)-zp(j,k-1,i-1)
     2   -zp(j+1,k-1,i-1))/dz
         g11im=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12im=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13im=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22im=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23im=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33im=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=0.25e0*(xp(j,k,i+1)+xp(j,k+1,i+1)+xp(j+1,k,i+1)
     1   +xp(j+1,k+1,i+1)-xp(j,k,i-1)-xp(j,k+1,i-1)-xp(j+1,k,i-1)
     2   -xp(j+1,k+1,i-1))/dx
         dypdx=0.25e0*(yp(j,k,i+1)+yp(j,k+1,i+1)+yp(j+1,k,i+1)
     1   +yp(j+1,k+1,i+1)-yp(j,k,i-1)-yp(j,k+1,i-1)-yp(j+1,k,i-1)
     2   -yp(j+1,k+1,i-1))/dx
         dzpdx=0.25e0*(zp(j,k,i+1)+zp(j,k+1,i+1)+zp(j+1,k,i+1)
     1   +zp(j+1,k+1,i+1)-zp(j,k,i-1)-zp(j,k+1,i-1)-zp(j+1,k,i-1)
     2   -zp(j+1,k+1,i-1))/dx
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k+1,i)-xp(j,k,i)
     1   -xp(j,k+1,i))/hy(j+1)
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k+1,i)-yp(j,k,i)
     1   -yp(j,k+1,i))/hy(j+1)
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k+1,i)-zp(j,k,i)
     1   -zp(j,k+1,i))/hy(j+1)
         dxpdz=0.5e0*(xp(j,k+1,i)+xp(j+1,k+1,i)-xp(j,k,i)
     1   -xp(j+1,k,i))/hz(k+1)
         dypdz=0.5e0*(yp(j,k+1,i)+yp(j+1,k+1,i)-yp(j,k,i)
     1   -yp(j+1,k,i))/hz(k+1)
         dzpdz=0.5e0*(zp(j,k+1,i)+zp(j+1,k+1,i)-zp(j,k,i)
     1   -zp(j+1,k,i))/hz(k+1)
         g11kp=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12kp=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13kp=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22kp=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23kp=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33kp=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dxpdx=0.25e0*(xp(j,k,i+1)+xp(j,k-1,i+1)+xp(j+1,k,i+1)
     1   +xp(j+1,k-1,i+1)-xp(j,k,i-1)-xp(j,k-1,i-1)-xp(j+1,k,i-1)
     2   -xp(j+1,k-1,i-1))/dx
         dypdx=0.25e0*(yp(j,k,i+1)+yp(j,k-1,i+1)+yp(j+1,k,i+1)
     1   +yp(j+1,k-1,i+1)-yp(j,k,i-1)-yp(j,k-1,i-1)-yp(j+1,k,i-1)
     2   -yp(j+1,k-1,i-1))/dx
         dzpdx=0.25e0*(zp(j,k,i+1)+zp(j,k-1,i+1)+zp(j+1,k,i+1)
     1   +zp(j+1,k-1,i+1)-zp(j,k,i-1)-zp(j,k-1,i-1)-zp(j+1,k,i-1)
     2   -zp(j+1,k-1,i-1))/dx
         dxpdy=0.5e0*(xp(j+1,k,i)+xp(j+1,k-1,i)-xp(j,k,i)
     1   -xp(j,k-1,i))/hy(j+1)
         dypdy=0.5e0*(yp(j+1,k,i)+yp(j+1,k-1,i)-yp(j,k,i)
     1   -yp(j,k-1,i))/hy(j+1)
         dzpdy=0.5e0*(zp(j+1,k,i)+zp(j+1,k-1,i)-zp(j,k,i)
     1   -zp(j,k-1,i))/hy(j+1)
         dxpdz=0.5e0*(xp(j,k,i)+xp(j+1,k,i)-xp(j,k-1,i)
     1   -xp(j+1,k-1,i))/hz(k)
         dypdz=0.5e0*(yp(j,k,i)+yp(j+1,k,i)-yp(j,k-1,i)
     1   -yp(j+1,k-1,i))/hz(k)
         dzpdz=0.5e0*(zp(j,k,i)+zp(j+1,k,i)-zp(j,k-1,i)
     1   -zp(j+1,k-1,i))/hz(k)
         g11km=dxpdx*dxpdx+dypdx*dypdx+dzpdx*dzpdx
         g12km=dxpdx*dxpdy+dypdx*dypdy+dzpdx*dzpdy
         g13km=dxpdx*dxpdz+dypdx*dypdz+dzpdx*dzpdz
         g22km=dxpdy*dxpdy+dypdy*dypdy+dzpdy*dzpdy
         g23km=dxpdy*dxpdz+dypdy*dypdz+dzpdy*dzpdz
         g33km=dxpdz*dxpdz+dypdz*dypdz+dzpdz*dzpdz
c     
         dgdx(j,k,1,1,1)=(g11ip-g11im)*2.0e0/dx
         dgdx(j,k,1,2,1)=(g12ip-g12im)*2.0e0/dx
         dgdx(j,k,1,3,1)=(g13ip-g13im)*2.0e0/dx
         dgdx(j,k,2,2,1)=(g22ip-g22im)*2.0e0/dx
         dgdx(j,k,2,3,1)=(g23ip-g23im)*2.0e0/dx
         dgdx(j,k,3,3,1)=(g33ip-g33im)*2.0e0/dx
         dgdx(j,k,2,1,1)=dgdx(j,k,1,2,1)
         dgdx(j,k,3,1,1)=dgdx(j,k,1,3,1)
         dgdx(j,k,3,2,1)=dgdx(j,k,2,3,1)
         dgdx(j,k,1,1,2)=(gv11(j+1,k)-gv11(j,k))/hy(j+1)
         dgdx(j,k,1,2,2)=(gv12(j+1,k)-gv12(j,k))/hy(j+1)
         dgdx(j,k,1,3,2)=(gv13(j+1,k)-gv13(j,k))/hy(j+1)
         dgdx(j,k,2,2,2)=(gv22(j+1,k)-gv22(j,k))/hy(j+1)
         dgdx(j,k,2,3,2)=(gv23(j+1,k)-gv23(j,k))/hy(j+1)
         dgdx(j,k,3,3,2)=(gv33(j+1,k)-gv33(j,k))/hy(j+1)
         dgdx(j,k,2,1,2)=dgdx(j,k,1,2,2)
         dgdx(j,k,3,1,2)=dgdx(j,k,1,3,2)
         dgdx(j,k,3,2,2)=dgdx(j,k,2,3,2)
         dgdx(j,k,1,1,3)=(g11kp-g11km)*2.0e0/dz
         dgdx(j,k,1,2,3)=(g12kp-g12km)*2.0e0/dz
         dgdx(j,k,1,3,3)=(g13kp-g13km)*2.0e0/dz
         dgdx(j,k,2,2,3)=(g22kp-g22km)*2.0e0/dz
         dgdx(j,k,2,3,3)=(g23kp-g23km)*2.0e0/dz
         dgdx(j,k,3,3,3)=(g33kp-g33km)*2.0e0/dz
         dgdx(j,k,2,1,3)=dgdx(j,k,1,2,3)
         dgdx(j,k,3,1,3)=dgdx(j,k,1,3,3)
         dgdx(j,k,3,2,3)=dgdx(j,k,2,3,3)
 150  continue
c     
      do 180 k=2,nz-1
      do 180 j=1,ny-1
         chf(j,k,2,1,1)=0.25e0*(gn12(j,k,i)+gn12(j+1,k,i))*
     1   (dgdx(j,k,1,1,1)+dgdx(j,k,1,1,1)-dgdx(j,k,1,1,1))
     2   +0.25e0*(gn22(j,k,i)+gn22(j+1,k,i))*
     3   (dgdx(j,k,1,2,1)+dgdx(j,k,2,1,1)-dgdx(j,k,1,1,2))
     4   +0.25e0*(gn23(j,k,i)+gn23(j+1,k,i))*
     5   (dgdx(j,k,1,3,1)+dgdx(j,k,3,1,1)-dgdx(j,k,1,1,3))
         chf(j,k,2,1,2)=0.25e0*(gn12(j,k,i)+gn12(j+1,k,i))*
     1   (dgdx(j,k,1,1,2)+dgdx(j,k,1,2,1)-dgdx(j,k,1,2,1))
     2   +0.25e0*(gn22(j,k,i)+gn22(j+1,k,i))*
     3   (dgdx(j,k,1,2,2)+dgdx(j,k,2,2,1)-dgdx(j,k,1,2,2))
     4   +0.25e0*(gn23(j,k,i)+gn23(j+1,k,i))*
     5   (dgdx(j,k,1,3,2)+dgdx(j,k,3,2,1)-dgdx(j,k,1,2,3))
         chf(j,k,2,1,3)=0.25e0*(gn12(j,k,i)+gn12(j+1,k,i))*
     1   (dgdx(j,k,1,1,3)+dgdx(j,k,1,3,1)-dgdx(j,k,1,3,1))
     2   +0.25e0*(gn22(j,k,i)+gn22(j+1,k,i))*
     3   (dgdx(j,k,1,2,3)+dgdx(j,k,2,3,1)-dgdx(j,k,1,3,2))
     4   +0.25e0*(gn23(j,k,i)+gn23(j+1,k,i))*
     5   (dgdx(j,k,1,3,3)+dgdx(j,k,3,3,1)-dgdx(j,k,1,3,3))
         chf(j,k,2,2,2)=0.25e0*(gn12(j,k,i)+gn12(j+1,k,i))*
     1   (dgdx(j,k,2,1,2)+dgdx(j,k,1,2,2)-dgdx(j,k,2,2,1))
     2   +0.25e0*(gn22(j,k,i)+gn22(j+1,k,i))*
     3   (dgdx(j,k,2,2,2)+dgdx(j,k,2,2,2)-dgdx(j,k,2,2,2))
     4   +0.25e0*(gn23(j,k,i)+gn23(j+1,k,i))*
     5   (dgdx(j,k,2,3,2)+dgdx(j,k,3,2,2)-dgdx(j,k,2,2,3))
         chf(j,k,2,2,3)=0.25e0*(gn12(j,k,i)+gn12(j+1,k,i))*
     1   (dgdx(j,k,2,1,3)+dgdx(j,k,1,3,2)-dgdx(j,k,2,3,1))
     2   +0.25e0*(gn22(j,k,i)+gn22(j+1,k,i))*
     3   (dgdx(j,k,2,2,3)+dgdx(j,k,2,3,2)-dgdx(j,k,2,3,2))
     4   +0.25e0*(gn23(j,k,i)+gn23(j+1,k,i))*
     5   (dgdx(j,k,2,3,3)+dgdx(j,k,3,3,2)-dgdx(j,k,2,3,3))
         chf(j,k,2,3,3)=0.25e0*(gn12(j,k,i)+gn12(j+1,k,i))*
     1   (dgdx(j,k,3,1,3)+dgdx(j,k,1,3,3)-dgdx(j,k,3,3,1))
     2   +0.25e0*(gn22(j,k,i)+gn22(j+1,k,i))*
     3   (dgdx(j,k,3,2,3)+dgdx(j,k,2,3,3)-dgdx(j,k,3,3,2))
     4   +0.25e0*(gn23(j,k,i)+gn23(j+1,k,i))*
     5   (dgdx(j,k,3,3,3)+dgdx(j,k,3,3,3)-dgdx(j,k,3,3,3))
         chf(j,k,2,2,1)=chf(j,k,2,1,2)
         chf(j,k,2,3,1)=chf(j,k,2,1,3)
         chf(j,k,2,3,2)=chf(j,k,2,2,3)
 180  continue
      return
      end



=END=