These changes should allow IDLinearWaves to compile...

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDLinearWaves/trunk@9 5c0f84ea-6048-4d6e-bfa4-55cd5f2e0dd7

4 files changed, 149 insertions, 145 deletions

diff --git a/param.ccl b/param.ccl
index 3acb784..5a8b69f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -3,8 +3,8 @@
 friend:einstein
 
 EXTENDS KEYWORD initial_data ""
-"teukwaves" : "linear waves initialdata: Teukolsky waves"
-"planewave" : "linear waves initialdata: plane waves"
+"teukwaves" :: "linear waves initialdata: Teukolsky waves"
+"planewave" :: "linear waves initialdata: plane waves"
 }
 
 private:
@@ -66,6 +66,6 @@ KEYWORD wavesgoing "in and outgoing waves..."
 
 KEYWORD teuk_no_vee "Initialize Teuk. waves with V=0?"
 {
-  "no" : " Bona Masso setting"
-  "yes": " Bona Masso setting"
+  "no" :: "Bona Masso setting"
+  "yes":: "Bona Masso setting"
 } "no"
diff --git a/schedule.ccl b/schedule.ccl
index b2f0d5d..f7b878f 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,6 +1,6 @@
 # Schedule definitions for thorn IDLinearWaves
 
-if (CCTK_Equals(initaldata,"linearwaves")) 
+if (CCTK_Equals(initial_data,"linearwaves")) 
 {
   schedule LinearWaves at CCTK_INITIAL
   {
diff --git a/src/planewaves.F b/src/planewaves.F
index 17b8431..683de4c 100644
--- a/src/planewaves.F
+++ b/src/planewaves.F
@@ -46,7 +46,7 @@
       
       INTEGER iin,iout
       REAL pi
-      REAL amplitude,ra,the,phi
+      REAL ra,the,phi
       REAL wave,wavep
       REAL kx,ky,kz,w
       REAL kxp,kyp,kzp,wp
@@ -54,11 +54,13 @@
 c     local arrays (scalars! Man, what a sweat!)
       REAL plus,minus,plusp,minusp,ain,aout
       
-      INTEGERE i,j,k
+      INTEGER i,j,k
+
+      INTEGER CCTK_Equals
 
       pi = 3.14159265358979d0
 
-      c     Wave characteristics
+c     Wave characteristics
 
 c     amplitude is declared by DECLARE_PARAMETER, for all other parameters
 c     we stick to the naming in the cactus3.2.0 version of cactus
@@ -106,7 +108,7 @@ c     precalc
       kzp = 2*pi*cos(the)/wavep
       wp = (kxp*kxp+kyp*kyp+kzp*kzp)**0.5
 
-      c  *************** plane waves ********************
+c  *************** plane waves ********************
       do k=1,sh(3)
          do j=1,sh(2)
             do i=1,sh(1)		 
@@ -141,7 +143,7 @@ c     the metric functions
 c     At this pint, the CAC3.2.0 version has some code for the dxgxx,....
 c     defined, how does CAC4.0 handles this without  #ifdef THORN_BONAMASSO ?
         
-        c  and finally the extrinsic curvatures
+c  and finally the extrinsic curvatures
 c       Joan: I removed the division by alp as the initialization
 c       of the lapse may take place later. This data is consistent
 c       for initial lapse equal one.
diff --git a/src/teukwaves.F b/src/teukwaves.F
index c7db62a..a2f1ca4 100644
--- a/src/teukwaves.F
+++ b/src/teukwaves.F
@@ -30,23 +30,22 @@
 
       subroutine teukwaves(CCTK_FARGUMENTS)
       implicit none
-
+      
       DECLARE_CCTK_FARGUMENTS
       DECLARE_PARAMETERS
-
+      
       
       REAL amp,m,ra,pi
       REAL wave,wave2,wave3,wave4,wave5,wave6,wave7,wave8
-      REAL wavepulse
       INTEGER ipacket,iparity
       INTEGER ingoing,outgoing
       REAL kappa
 c     point values of x.y z,r
       REAL xp,yp,zp,rp
-
+      
 c     spherical metric
-      REAL grr,grt,grp,gtt,gtp,gpp
-      REAL hrr,hrt,hrp,htt,htp,hpp
+      REAL teuk_grr,teuk_grt,teuk_grp,teuk_gtt,teuk_gtp,teuk_gpp
+      REAL teuk_hrr,teuk_hrt,teuk_hrp,teuk_htt,teuk_htp,teuk_hpp
       
 c     special coefficeints needed for teukolsky waves
 c     All these were 3d arrays in old versions of Newage,G,etc..
@@ -61,14 +60,17 @@ c     from old spheretocart
       REAL drx,dry,drz,dtx,dty,dtz,dpx,dpy,dpz
       
       INTEGER i,j,k
-
+      
+      INTEGER CCTK_Equals
+      
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
       pi    = 3.14159265358979
       kappa = sqrt(143.0d0/pi)/12288.0d0
-
+      
       if (wavecenter.ne.0.) then
-         CCTK_Warn(0,"Teukwaves: wavecenter must be zero for time symmetry")
+         call CCTK_Warn(0,"Teukwaves: wavecenter must be 
+     &   zero for time symmetry")
       endif
  
 c     CACTUS 4.0: all REAL/INTEGER parameter are available as variables and 
@@ -76,24 +78,24 @@ c     need not be derived from the database. Nevertheless I will assign
 c     the parameter variables to the variables of CAC3.0, to allow for 
 c     easy portation (wavepulse stays wavepulse btw). (GERD)
       
-      am  = amplitude
+      amp  = amplitude
       m   = mvalue
       wave= wavelength
-
+      
 c     to be consistent with Evans & Abrahams, rescale the amplitude
 c     by kappa, because they use that in the packet and
 c     4/3 because they use I not Teuk. F  
-      if (CCTK_Equals(packet,'eppley')) then
+      if (CCTK_Equals(packet,'eppley') == 1) then
          ipacket = 1
          write(*,*)'Teuk.Packet= eppley' 
-      elseif (CCTK_Equals(packet,'evans')) then
+      elseif (CCTK_Equals(packet,'evans') == 1) then
          ipacket = 2
          write(*,*)'Teuk.Packet= evans'
 c     to be consistent with Evans & Abrahams, rescale the amplitude
 c     by kappa, because they use that in the packet and
 c     4/3 because they use I not Teuk. F
          amp = amp*(4.0d0/3.0d0)*kappa
-      elseif (CCTK_Equals(packet,'square')) then
+      elseif (CCTK_Equals(packet,'square') == 1) then
          ipacket = 3
          write(*,*)'Teuk.Packet= square'
       else
@@ -107,25 +109,25 @@ c     4/3 because they use I not Teuk. F
       write(*,*)'wavecenter = ',ra
       write(*,*)'wavelength = ',wave
       
-      if (CCTK_Equals(parity,'even')) then    
+      if (CCTK_Equals(parity,'even') == 1) then    
          iparity = 0
-      else (CONTAINS(parity,'odd')) then 
+      elseif (CCTK_Equals(parity,'odd') == 1) then 
          iparity = 1
       endif
-
-       if(CCTK_Equals(wavegoing,'in')) then
-          ingoing = 1
-       endif
-
-      if(CCTK_Equals(wavegoing,'out')) then
+      
+      if(CCTK_Equals(wavesgoing,'in') == 1) then
+         ingoing = 1
+      endif
+      
+      if(CCTK_Equals(wavesgoing,'out') == 1) then
          outgoing = 1
       endif
-
-      if(CCTK_Equals(wavegoing,'both')) then
+      
+      if(CCTK_Equals(wavesgoing,'both') == 1) then
          ingoing  = 1
          outgoing = 1
       endif
-
+      
       if (ipacket.eq.1) then
          if (ingoing.eq.0.or.outgoing.eq.0) then
             write(*,*)'WARNING: Epply packet is only non-singular at'
@@ -183,9 +185,9 @@ c     to keep the analytic solution valid for the eppley packet:
                      teuk_tn = (time-rp-ra)
                      tn2 = teuk_tn**2
                   endif
-
+                  
                   if (ingoing.eq.1) then
-
+                     
                      tmp = exp(-tp2)
                      fp = amp * teuk_tp*tmp
                      fpa = amp * (1 - 2*tp2)*tmp
@@ -194,23 +196,23 @@ c     to keep the analytic solution valid for the eppley packet:
                      fpd = amp * teuk_tp*(60 - 80*tp2 + 16*tp2*tp2)*tmp
                      fpe = amp * (
      &                    (60 - 360*tp2 + 240*tp2*tp2 - 32*tp2*tp2*tp2)*tmp )
-
+                     
                      if (outgoing.eq.0.and.(teuk_tp+wavepulse).le.0) 
      &                    then
-
+                        
                         fp  = 0.0d0
                         fpa = 0.0d0
                         fpb = 0.0d0
                         fpc = 0.0d0
                         fpd = 0.0d0
                         fpe = 0.0d0
-
+                        
                      endif
-
+                     
                   endif
-
+                  
                   if (outgoing.eq.1) then
-
+                     
                      tmp = exp(-tn2)
                      fn = amp * teuk_tn*tmp
                      fna = amp * (1 - 2*tn2)*tmp
@@ -219,36 +221,36 @@ c     to keep the analytic solution valid for the eppley packet:
                      fnd = amp * teuk_tn*(60 - 80*tn2 + 16*tn2*tn2)*tmp
                      fne = amp * (
      &                    (60 - 360*tn2 + 240*tn2*tn2 - 32*tn2*tn2*tn2)*tmp )
-
+                     
                   endif
                endif
-
-      
+               
+               
 c     evans package  --> w^4(1-x^2/w^2)^6
                if(ipacket.eq.2) then
-
+                  
                   if (ingoing.eq.1) then
-
+                     
                      tmp = 1.0d0 - tp2/wave2
-
+                     
                      if (abs(teuk_tp).lt.wave) then
-
+                        
                         fp  =  amp*wave5*tmp**6
                         fpa =  -12.0d0*amp*wave3*teuk_tp*tmp**5
                         fpb =  -12.0d0*amp*wave3*(
      &                       tmp**5 - 10.0d0*tp2/wave2*tmp**4 )
                         fpc =  120.0d0*amp*wave3*(
-     &                        3.0d0*teuk_tp/wave2*tmp**4 
+     &                       3.0d0*teuk_tp/wave2*tmp**4 
      &                       -8.0d0*teuk_tp*tp2/wave4*tmp**3 )
                         fpd = 360.0d0*amp*wave3*(
-     &                        1.0d0/wave2*tmp**4 
+     &                       1.0d0/wave2*tmp**4 
      &                       -16.0d0*tp2/wave4*tmp**3
      &                       +16.0d0*tp2*tp2/wave6*tmp**2 )
                         fpe = 2880.0d0*amp*wave3*(
      &                       -5.0d0*teuk_tp/wave4*tmp**3
      &                       +20.0d0*teuk_tp*tp2/wave6*tmp**2
      &                       -8.0d0*teuk_tp*tp2*tp2/wave8*tmp )
-
+                        
                      else
                         fp=0.
                         fpa=0.
@@ -257,31 +259,31 @@ c     evans package  --> w^4(1-x^2/w^2)^6
                         fpd=0.
                         fpe=0.
                      endif
-
+                     
                   endif
-
+                  
                   if (outgoing.eq.1) then
-
+                     
                      tmp = 1.0d0 - tn2/wave2
-         
+                     
                      if (abs(teuk_tn).lt.wave) then
-
+                        
                         fn  =  amp*wave5*tmp**6
                         fna =  -12.0d0*amp*wave3*teuk_tn*tmp**5
                         fnb =  -12.0d0*amp*wave3*(
      &                       tmp**5 - 10.0d0*tn2/wave2*tmp**4 )
                         fnc =  120.0d0*amp*wave3*(
-     &                        3.0d0*teuk_tn/wave2*tmp**4 
+     &                       3.0d0*teuk_tn/wave2*tmp**4 
      &                       -8.0d0*teuk_tn*tn2/wave4*tmp**3 )
                         fnd = 360.0d0*amp*wave3*(
-     &                        1.0d0/wave2*tmp**4 
+     &                       1.0d0/wave2*tmp**4 
      &                       -16.0d0*tn2/wave4*tmp**3
      &                       +16.0d0*tn2*tn2/wave6*tmp**2 )
                         fne = 2880.0d0*amp*wave3*(
      &                       -5.0d0*teuk_tn/wave4*tmp**3
      &                       +20.0d0*teuk_tn*tn2/wave6*tmp**2
      &                       -8.0d0*teuk_tn*tn2*tn2/wave8*tmp )
-
+                        
                      else
                         fn=0.
                         fna=0.
@@ -290,17 +292,17 @@ c     evans package  --> w^4(1-x^2/w^2)^6
                         fnd=0.
                         fne=0.
                      endif
-
+                     
                   endif
-
+                  
                endif
       
 c     square package --> (1-x^2/w^2)^2
                if(ipacket.eq.3) then 
-
+                  
                   write(*,*)'WARNING: need to calculate fpe,fne'
                   write(*,*)'WARNING: need conditionals for in out waves'
-
+                  
                   if (abs(teuk_tp).lt.wave) then
                      fp = amp * (1-(tp2/wave2))**2
                      fpa = amp * (-4*teuk_tp)*(1-tp2/wave2)/wave2
@@ -329,7 +331,7 @@ c     square package --> (1-x^2/w^2)^2
                      fnd=0.
                   endif
                endif
-      
+               
 c     
 c     coefficients
 c     
@@ -347,7 +349,7 @@ c     of the above coefficients.  That is because the above coefficients
 c     assume an ingoing and outgoing wave packet of the form:
 c     wave = ingoing - outgoing.  For a single wave, we want ingoing
 c     and outgoing to have the same sign.
-
+               
                if(ingoing.eq.0.and.outgoing.eq.1) then
                   ca=-ca
                   cb=-cb
@@ -438,23 +440,23 @@ c     mvalue
                
 c     parity
                if (iparity.eq.0) then    
-                  grr = 1 + ca*frr
-                  grt = cb*frt*rp
-                  grp = cb*frp*(xp**2+yp**2)**0.5
-                  gtt = rp**2*(1 + cc*ftt1 + ca*ftt2)
-                  gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5
-                  gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2)
+                  teuk_grr = 1 + ca*frr
+                  teuk_grt = cb*frt*rp
+                  teuk_grp = cb*frp*(xp**2+yp**2)**0.5
+                  teuk_gtt = rp**2*(1 + cc*ftt1 + ca*ftt2)
+                  teuk_gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5
+                  teuk_gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2)
                else
-                  grr = 1.
-                  grt = ck*drt*rp
-                  grp = ck*drp*rp*sint
-                  gtt = (1+cl*dtt)*rp**2
-                  gtp = cl*dtp*rp**2*sint
-                  gpp = (1+cl*dpp)*rp**2*sint**2
+                  teuk_grr = 1.
+                  teuk_grt = ck*drt*rp
+                  teuk_grp = ck*drp*rp*sint
+                  teuk_gtt = (1+cl*dtt)*rp**2
+                  teuk_gtp = cl*dtp*rp**2*sint
+                  teuk_gpp = (1+cl*dpp)*rp**2*sint**2
                endif
-
+               
 c     Finally, convert the spherical components to cartesian
-
+               
 c     define derivatives drx = (dr/dx)               
                drx = xp/rp
                dry = yp/rp
@@ -467,31 +469,34 @@ c     define derivatives drx = (dr/dx)
                dpx = -yp/(xp**2+yp**2)
                dpy = xp/(xp**2+yp**2)
                dpz = 0.
-
+               
 c     cartesian g
-               gxx(i,j,k) = drx**2*grr+2.*drx*dtx*grt
-     &              +2.*drx*dpx*grp+dtx**2*gtt
-     &              +2.*dtx*dpx*gtp+dpx**2*gpp
+               gxx(i,j,k) = drx**2*teuk_grr+2.*drx*dtx*teuk_grt
+     &              +2.*drx*dpx*teuk_grp+dtx**2*teuk_gtt
+     &              +2.*dtx*dpx*teuk_gtp+dpx**2*teuk_gpp
                
-               gyy(i,j,k) = dry**2*grr+2.*dry*dty*grt
-     &              +2.*dry*dpy*grp+dty**2*gtt
-     &              +2.*dty*dpy*gtp+dpy**2*gpp
+               gyy(i,j,k) = dry**2*teuk_grr+2.*dry*dty*teuk_grt
+     &              +2.*dry*dpy*teuk_grp+dty**2*teuk_gtt
+     &              +2.*dty*dpy*teuk_gtp+dpy**2*teuk_gpp
  
-               gzz(i,j,k) = drz**2*grr+2.*drz*dtz*grt
-     &              +2.*drz*dpz*grp+dtz**2*gtt
-     &              +2.*dtz*dpz*gtp+dpz**2*gpp
+               gzz(i,j,k) = drz**2*teuk_grr+2.*drz*dtz*teuk_grt
+     &              +2.*drz*dpz*teuk_grp+dtz**2*teuk_gtt
+     &              +2.*dtz*dpz*teuk_gtp+dpz**2*teuk_gpp
                
-               gxy(i,j,k) = drx*dry*grr+(drx*dty+dtx*dry)*grt
-     &              +(drx*dpy+dpx*dry)*grp
-     &              +dtx*dty*gtt+(dtx*dpy+dpx*dty)*gtp+dpx*dpy*gpp
+               gxy(i,j,k) = drx*dry*teuk_grr+(drx*dty+dtx*dry)*teuk_grt
+     &              +(drx*dpy+dpx*dry)*teuk_grp
+     &              +dtx*dty*teuk_gtt+(dtx*dpy+dpx*dty)*teuk_gtp+
+     &              dpx*dpy*teuk_gpp
                
-               gxz(i,j,k) = drx*drz*grr+(drx*dtz+dtx*drz)*grt
-     &              +(drx*dpz+dpx*drz)*grp
-     &              +dtx*dtz*gtt+(dtx*dpz+dpx*dtz)*gtp+dpx*dpz*gpp
+               gxz(i,j,k) = drx*drz*teuk_grr+(drx*dtz+dtx*drz)*teuk_grt
+     &              +(drx*dpz+dpx*drz)*teuk_grp
+     &              +dtx*dtz*teuk_gtt+(dtx*dpz+dpx*dtz)*teuk_gtp+
+     &              dpx*dpz*teuk_gpp
                
-               gyz(i,j,k) = dry*drz*grr+(dry*dtz+dty*drz)*grt
-     &              +(dry*dpz+dpy*drz)*grp
-     &              +dty*dtz*gtt+(dty*dpz+dpy*dtz)*gtp+dpy*dpz*gpp
+               gyz(i,j,k) = dry*drz*teuk_grr+(dry*dtz+dty*drz)*teuk_grt
+     &              +(dry*dpz+dpy*drz)*teuk_grp
+     &              +dty*dtz*teuk_gtt+(dty*dpz+dpy*dtz)*teuk_gtp+
+     &              dpy*dpz*teuk_gpp
                 
 
 
@@ -506,64 +511,61 @@ c     for the ext. curv. well need...
                ckdot = (fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3+3*(fpa-fna)/rp**4
                cldot = -(fpd+fnd)/rp+2*(fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3+
      &              3*(fpa-fna)/rp**4
-
+               
 c     parity
                if (iparity.eq.0) then    
-                  hrr = -0.5*(cadot*frr)
-                  hrt = -0.5*(cbdot*frt*rp)
-                  hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5)
-                  htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2))
-                  htp = -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5)
-                  hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp2))
+                  teuk_hrr = -0.5*(cadot*frr)
+                  teuk_hrt = -0.5*(cbdot*frt*rp)
+                  teuk_hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5)
+                  teuk_htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2))
+                  teuk_htp = 
+     &                 -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5)
+                  teuk_hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp2))
                else
-                  hrr = 0.
-                  hrt = -0.5*(ckdot*drt*rp)
-                  hrp = -0.5*(ckdot*drp*rp*sint)
-                  htt = -0.5*((cldot*dtt)*rp**2)
-                  htp = -0.5*(cldot*dtp*rp**2*sint)
-                  hpp = -0.5*((cldot*dpp)*rp**2*sint**2)
+                  teuk_hrr = 0.
+                  teuk_hrt = -0.5*(ckdot*drt*rp)
+                  teuk_hrp = -0.5*(ckdot*drp*rp*sint)
+                  teuk_htt = -0.5*((cldot*dtt)*rp**2)
+                  teuk_htp = -0.5*(cldot*dtp*rp**2*sint)
+                  teuk_hpp = -0.5*((cldot*dpp)*rp**2*sint**2)
                endif
                
 c     time symmetry
-               hxx(i,j,k) = drx**2*hrr+2.*drx*dtx*hrt
-     &              +2.*drx*dpx*hrp+dtx**2*htt
-     &              +2.*dtx*dpx*htp+dpx**2*hpp
-
-               hyy(i,j,k) = dry**2*hrr+2.*dry*dty*hrt
-     &              +2.*dry*dpy*hrp+dty**2*htt
-     &              +2.*dty*dpy*htp+dpy**2*hpp
+               hxx(i,j,k) = drx**2*teuk_hrr+2.*drx*dtx*teuk_hrt
+     &              +2.*drx*dpx*teuk_hrp+dtx**2*teuk_htt
+     &              +2.*dtx*dpx*teuk_htp+dpx**2*teuk_hpp
                
-               hzz(i,j,k) = drz**2*hrr+2.*drz*dtz*hrt
-     &              +2.*drz*dpz*hrp+dtz**2*htt
-     &              +2.*dtz*dpz*htp+dpz**2*hpp
+               hyy(i,j,k) = dry**2*teuk_hrr+2.*dry*dty*teuk_hrt
+     &              +2.*dry*dpy*teuk_hrp+dty**2*teuk_htt
+     &              +2.*dty*dpy*teuk_htp+dpy**2*teuk_hpp
                
-               hxy(i,j,k) = drx*dry*hrr+(drx*dty+dtx*dry)*hrt
-     &              +(drx*dpy+dpx*dry)*hrp
-     &              +dtx*dty*htt+(dtx*dpy+dpx*dty)*htp+dpx*dpy*hpp
+               hzz(i,j,k) = drz**2*teuk_hrr+2.*drz*dtz*teuk_hrt
+     &              +2.*drz*dpz*teuk_hrp+dtz**2*teuk_htt
+     &              +2.*dtz*dpz*teuk_htp+dpz**2*teuk_hpp
+               
+               hxy(i,j,k) = drx*dry*teuk_hrr+(drx*dty+dtx*dry)*teuk_hrt
+     &              +(drx*dpy+dpx*dry)*teuk_hrp
+     &              +dtx*dty*teuk_htt+(dtx*dpy+dpx*dty)*teuk_htp+
+     &              dpx*dpy*teuk_hpp
+               
+               hxz(i,j,k) = drx*drz*teuk_hrr+(drx*dtz+dtx*drz)*teuk_hrt
+     &              +(drx*dpz+dpx*drz)*teuk_hrp
+     &              +dtx*dtz*teuk_htt+(dtx*dpz+dpx*dtz)*teuk_htp+
+     &              dpx*dpz*teuk_hpp
+               
+               hyz(i,j,k) = dry*drz*teuk_hrr+(dry*dtz+dty*drz)*teuk_hrt
+     &              +(dry*dpz+dpy*drz)*teuk_hrp
+     &              +dty*dtz*teuk_htt+(dty*dpz+dpy*dtz)*teuk_htp+
+     &              dpy*dpz*teuk_hpp
                
-               hxz(i,j,k) = drx*drz*hrr+(drx*dtz+dtx*drz)*hrt
-     &              +(drx*dpz+dpx*drz)*hrp
-     &              +dtx*dtz*htt+(dtx*dpz+dpx*dtz)*htp+dpx*dpz*hpp
                
-               hyz(i,j,k) = dry*drz*hrr+(dry*dtz+dty*drz)*hrt
-     &              +(dry*dpz+dpy*drz)*hrp
-     &              +dty*dtz*htt+(dty*dpz+dpy*dtz)*htp+dpy*dpz*hpp
-                
-
                
-c$$$               hxx(i,j,k) = 0.
-c$$$               hxy(i,j,k) = 0.
-c$$$               hxz(i,j,k) = 0.
-c$$$               hyy(i,j,k) = 0.
-c$$$               hyz(i,j,k) = 0.
-c$$$               hzz(i,j,k) = 0.
-
-
             enddo
          enddo
       enddo
-
- 
+      
       return
-      end
+      end 
+      
+