      module types
      integer ,  parameter:: I4B=SELECTED_INT_KIND (9) 
      integer ,  parameter:: I2B=SELECTED_INT_KIND (4) 
      integer ,  parameter:: I1B=SELECTED_INT_KIND (2) 
      integer ,  parameter:: SP=KIND (1.0D0) 
      end module


      module hpfprocessors
!HPF$      processors PROC(4)
      end module


      module diff_interface
      interface 
         subroutine diff_x(nx, ny, nz, h, f, df)
         use types
         use hpfprocessors
         integer (I4B):: nx,ny,nz
         real (SP):: h
         real (SP),  dimension(:,:,:),intent(in):: f 
         real (SP),  dimension(:,:,:),intent(out):: df 
!HPF$      distribute * (*,*,block) onto * PROC :: f
!HPF$      align with * f :: df
         end subroutine

         subroutine diff_y(nx, ny, nz, h, f, df)
         use types
         use hpfprocessors
         integer (I4B):: nx,ny,nz
         real (SP):: h
         real (SP),  dimension(:,:,:),intent(in):: f 
         real (SP),  dimension(:,:,:),intent(out):: df 
!HPF$      distribute * (*,*,block) onto * PROC :: f
!HPF$      align with * f :: df
         end subroutine

         subroutine diff_z(nx, ny, nz, h, f, df)
         use types
         use hpfprocessors
         integer (I4B):: nx,ny,nz
         real (SP):: h
         real (SP),  dimension(:,:,:),intent(in):: f 
         real (SP),  dimension(:,:,:),intent(out):: df 
!HPF$      distribute * (*,*,block) onto * PROC :: f
!HPF$      align with * f :: df
         end subroutine

      end interface
      end module


      module evol_interface
      interface 
         subroutine tensors(g11, g12, g13, g22, g23, g33, gup11, gup12, &
     &gup13, gup22, gup23, gup33, c111, c112, c113, c122, c123, c133, c2&
     &11, c212, c213, c222, c223, c233, c311, c312, c313, c322, c323, c3&
     &33, r11, r12, r13, r22, r23, r33, nx, ny, nz, h)
         use hpfprocessors
         use types
         real (SP),  dimension(:,:,:),intent(in):: g11,g12,g13,g22,g23,g&
     &33 
         real (SP),  dimension(:,:,:),intent(out):: gup11,gup12,gup13,gu&
     &p22,gup23,gup33 
         real (SP),  dimension(:,:,:),intent(out):: c111,c112,c113,c122,&
     &c123,c133 
         real (SP),  dimension(:,:,:),intent(out):: c211,c212,c213,c222,&
     &c223,c233 
         real (SP),  dimension(:,:,:),intent(out):: c311,c312,c313,c322,&
     &c323,c333 
         real (SP),  dimension(:,:,:),intent(out):: r11,r12,r13,r22,r23,&
     &r33 
!HPF$      distribute g11 * (*,*,block) onto * PROC
!HPF$      align with * g11 :: g12,g13,g22,g23,g33
!HPF$      align with * g11 :: gup11,gup12,gup13,gup22,gup23,gup33
!HPF$      align with * g11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with * g11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with * g11 :: c311,c312,c313,c322,c323,c333
!HPF$      align with * g11 :: r11,r12,r13,r22,r23,r33
         integer (I4B):: nx,ny,nz
         real (SP):: h
         end subroutine

         subroutine evolve_g(g11, g12, g13, g22, g23, g33, g11m1, g12m1,&
     & g13m1, g22m1, g23m1, g33m1, c111, c112, c113, c122, c123, c133, c&
     &211, c212, c213, c222, c223, c233, c311, c312, c313, c322, c323, c&
     &333, k11, k12, k13, k21, k22, k23, k31, k32, k33, alpha, beta1, be&
     &ta2, beta3, h, dt, nx, ny, nz)
         use types
         use hpfprocessors
         use diff_interface
         real (SP),  dimension(:,:,:),intent(in):: g11,g12,g13,g22,g23,g&
     &33 
         real (SP),  dimension(:,:,:),intent(out):: g11m1,g12m1,g13m1,g2&
     &2m1,g23m1,g33m1 
         real (SP),  dimension(:,:,:),intent(in):: c111,c112,c113,c122,c&
     &123,c133 
         real (SP),  dimension(:,:,:),intent(in):: c211,c212,c213,c222,c&
     &223,c233 
         real (SP),  dimension(:,:,:),intent(in):: c311,c312,c313,c322,c&
     &323,c333 
         real (SP),  dimension(:,:,:),intent(in):: k11,k12,k13,k21,k22,k&
     &23,k31,k32,k33 
         real (SP),  dimension(:,:,:),intent(in):: alpha,beta1,beta2,bet&
     &a3 
         integer (I4B):: nx,ny,nz
         real (SP):: h,dt
!HPF$      distribute g11 * (*,*,block) onto * PROC
!HPF$      align with * g11 :: g12,g13,g22,g23,g33
!HPF$      align with * g11 :: g11m1,g12m1,g13m1,g22m1,g23m1,g33m1
!HPF$      align with * g11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with * g11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with * g11 :: c311,c312,c313,c322,c323,c333
!HPF$      align with * g11 :: k11,k12,k13,k21,k22,k23,k31,k32,k33
!HPF$      align with * g11 :: alpha,beta1,beta2,beta3
         end subroutine

         subroutine evolve_k(gup11, gup12, gup13, gup22, gup23, gup33, k&
     &11, k12, k13, k21, k22, k23, k31, k32, k33, k11m1, k12m1, k13m1, k&
     &21m1, k22m1, k23m1, k31m1, k32m1, k33m1, c111, c112, c113, c122, c&
     &123, c133, c211, c212, c213, c222, c223, c233, c311, c312, c313, c&
     &322, c323, c333, r11, r12, r13, r22, r23, r33, alpha, beta1, beta2&
     &, beta3, nx, ny, nz, h, dt)
         use types
         use hpfprocessors
         real (SP),  dimension(:,:,:),intent(in):: gup11,gup12,gup13,gup&
     &22,gup23,gup33 
         real (SP),  dimension(:,:,:),intent(in):: k11,k12,k13,k21,k22,k&
     &23,k31,k32,k33 
         real (SP),  dimension(:,:,:),intent(inout):: k11m1,k12m1,k13m1,&
     &k21m1,k22m1,k23m1,k31m1,k32m1,k33m1 
         real (SP),  dimension(:,:,:),intent(in):: c111,c112,c113,c122,c&
     &123,c133 
         real (SP),  dimension(:,:,:),intent(in):: c211,c212,c213,c222,c&
     &223,c233 
         real (SP),  dimension(:,:,:),intent(in):: c311,c312,c313,c322,c&
     &323,c333 
         real (SP),  dimension(:,:,:),intent(in):: r11,r12,r13,r22,r23,r&
     &33 
         real (SP),  dimension(:,:,:),intent(in):: alpha,beta1,beta2,bet&
     &a3 
!HPF$      distribute * (*,*,block) onto * PROC :: k11
!HPF$      align with * k11 :: k12,k13,k21,k22,k23,k31,k32,k33
!HPF$      align with * k11 :: k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1&
!HPF$&,k32m1,k33m1
!HPF$      align with * k11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with * k11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with * k11 :: c311,c312,c313,c322,c323,c333
!HPF$      align with * k11 :: r11,r12,r13,r22,r23,r33
!HPF$      align with * k11 :: gup11,gup12,gup13,gup22,gup23,gup33
!HPF$      align with * k11 :: alpha,beta1,beta2,beta3
         integer (I4B):: nx,ny,nz
         real (SP):: h,dt
         end subroutine

         subroutine swap(a, b)
         use types
         use hpfprocessors
         real (SP),  dimension(:,:,:),intent(inout):: a,b 
!HPF$      distribute * (*,*,block) onto * PROC :: a
!HPF$      align with * a :: b
         end subroutine

         subroutine constraint(nx, ny, nz, k11, k12, k13, k21, k22, k23,&
     & k31, k32, k33, c111, c112, c113, c122, c123, c133, c211, c212, c2&
     &13, c222, c223, c233, c311, c312, c313, c322, c323, c333, r11, r12&
     &, r13, r22, r23, r33, gup11, gup12, gup13, gup22, gup23, gup33, ha&
     &mnorm, mom1norm, mom2norm, mom3norm, h)
         use types
         use hpfprocessors
         use diff_interface
         integer (I4B):: nx,ny,nz
         real (SP),  dimension(:,:,:),intent(in):: k11,k12,k13 
         real (SP),  dimension(:,:,:),intent(in):: k21,k22,k23 
         real (SP),  dimension(:,:,:),intent(in):: k31,k32,k33 
         real (SP),  dimension(:,:,:),intent(in):: r11,r12,r13 
         real (SP),  dimension(:,:,:),intent(in):: r22,r23 
         real (SP),  dimension(:,:,:),intent(in):: r33 
         real (SP),  dimension(:,:,:),intent(in):: gup11,gup12,gup13 
         real (SP),  dimension(:,:,:),intent(in):: gup22,gup23 
         real (SP),  dimension(:,:,:),intent(in):: gup33 
         real (SP),  dimension(:,:,:),intent(in):: c111,c112,c113 
         real (SP),  dimension(:,:,:),intent(in):: c122,c123 
         real (SP),  dimension(:,:,:),intent(in):: c133 
         real (SP),  dimension(:,:,:),intent(in):: c211,c212,c213 
         real (SP),  dimension(:,:,:),intent(in):: c222,c223 
         real (SP),  dimension(:,:,:),intent(in):: c233 
         real (SP),  dimension(:,:,:),intent(in):: c311,c312,c313 
         real (SP),  dimension(:,:,:),intent(in):: c322,c323 
         real (SP),  dimension(:,:,:),intent(in):: c333 
!HPF$      distribute * (*,*,block) onto * PROC :: k11
!HPF$      align with * k11 :: k12,k13,k21,k22,k23,k31,k32,k33
!HPF$      align with * k11 :: r11,r12,r13,r22,r23,r33
!HPF$      align with * k11 :: gup11,gup12,gup13,gup22,gup23,gup33
!HPF$      align with * k11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with * k11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with * k11 :: c311,c312,c313,c322,c323,c333
         real (SP):: hamnorm,mom1norm,mom2norm,mom3norm,h
         end subroutine

      end interface
      end module


      subroutine initialize(g11, g12, g13, g22, g23, g33, g11m1, g12m1, &
     &g13m1, g22m1, g23m1, g33m1, k11, k12, k13, k21, k22, k23, k31, k32&
     &, k33, k11m1, k12m1, k13m1, k21m1, k22m1, k23m1, k31m1, k32m1, k33&
     &m1, alpha, beta1, beta2, beta3, nx, ny, nz, h, dt, time)
      use types
      use hpfprocessors
      use evol_interface
      real (SP),  dimension(:,:,:),intent(inout):: g11,g12,g13,g22,g23,g&
     &33 
      real (SP),  dimension(:,:,:),intent(inout):: g11m1,g12m1,g13m1,g22&
     &m1,g23m1,g33m1 
      real (SP),  dimension(:,:,:),intent(inout):: k11,k12,k13,k21,k22,k&
     &23,k31,k32,k33 
      real (SP),  dimension(:,:,:),intent(inout):: k11m1,k12m1,k13m1,k21&
     &m1,k22m1,k23m1,k31m1,k32m1,k33m1 
      real (SP),  dimension(:,:,:),intent(in):: alpha,beta1,beta2,beta3 

!HPF$      distribute g11 * (*,*,block) onto * PROC
!HPF$      align with * g11 :: g12,g13,g22,g23,g33
!HPF$      align with * g11 :: g11m1,g12m1,g13m1,g22m1,g23m1,g33m1
!HPF$      align with * g11 :: k11,k12,k13,k21,k22,k23,k31,k32,k33
!HPF$      align with * g11 :: k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1&
!HPF$&,k32m1,k33m1
!HPF$      align with * g11 :: alpha,beta1,beta2,beta3
      integer (I4B):: nx,ny,nz
      real (SP):: h,time,dt
      real (SP),  dimension(size (g11,1),size (g11,2),size (g11,3)):: gu&
     &p11,gup12,gup13,gup22,gup23,gup33,c111,c112,c113,c122,c123,c133,c2&
     &11,c212,c213,c222,c223,c233,c311,c312,c313,c322,c323,c333,r11,r12,&
     &r13,r22,r23,r33 
!HPF$      align with g11 :: gup11,gup12,gup13,gup22,gup23,gup33
!HPF$      align with g11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with g11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with g11 :: c311,c312,c313,c322,c323,c333
!HPF$      align with g11 :: r11,r12,r13,r22,r23,r33
      real (SP):: dtsave,timesave
      call tensors (g11,g12,g13,g22,g23,g33,gup11,gup12,gup13,gup22,gup2&
     &3,gup33,c111,c112,c113,c122,c123,c133,c211,c212,c213,c222,c223,c23&
     &3,c311,c312,c313,c322,c323,c333,r11,r12,r13,r22,r23,r33,nx,ny,nz,h&
     &)
      dtsave = dt
      timesave = time
      dt = (-(0.5d0*dt))
      dt = 0.5d0*dt
      g11m1 = g11
      g12m1 = g12
      g13m1 = g13
      g22m1 = g22
      g23m1 = g23
      g33m1 = g33
      k11m1 = k11
      k12m1 = k12
      k13m1 = k13
      k21m1 = k21
      k22m1 = k22
      k23m1 = k23
      k31m1 = k31
      k32m1 = k32
      k33m1 = k33
      call evolve_g (g11,g12,g13,g22,g23,g33,g11m1,g12m1,g13m1,g22m1,g23&
     &m1,g33m1,c111,c112,c113,c122,c123,c133,c211,c212,c213,c222,c223,c2&
     &33,c311,c312,c313,c322,c323,c333,k11,k12,k13,k21,k22,k23,k31,k32,k&
     &33,alpha,beta1,beta2,beta3,h,dt,nx,ny,nz)
      call evolve_k (gup11,gup12,gup13,gup22,gup23,gup33,k11,k12,k13,k21&
     &,k22,k23,k31,k32,k33,k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32&
     &m1,k33m1,c111,c112,c113,c122,c123,c133,c211,c212,c213,c222,c223,c2&
     &33,c311,c312,c313,c322,c323,c333,r11,r12,r13,r22,r23,r33,alpha,bet&
     &a1,beta2,beta3,nx,ny,nz,h,dt)
      call swap (g11,g11m1)
      call swap (g12,g12m1)
      call swap (g13,g13m1)
      call swap (g22,g22m1)
      call swap (g23,g23m1)
      call swap (g33,g33m1)
      call swap (k11,k11m1)
      call swap (k12,k12m1)
      call swap (k13,k13m1)
      call swap (k21,k21m1)
      call swap (k22,k22m1)
      call swap (k23,k23m1)
      call swap (k31,k31m1)
      call swap (k32,k32m1)
      call swap (k33,k33m1)
      dt = 2.0d0*dt
      time = time+dt
      call tensors (g11,g12,g13,g22,g23,g33,gup11,gup12,gup13,gup22,gup2&
     &3,gup33,c111,c112,c113,c122,c123,c133,c211,c212,c213,c222,c223,c23&
     &3,c311,c312,c313,c322,c323,c333,r11,r12,r13,r22,r23,r33,nx,ny,nz,h&
     &)
      call evolve_g (g11,g12,g13,g22,g23,g33,g11m1,g12m1,g13m1,g22m1,g23&
     &m1,g33m1,c111,c112,c113,c122,c123,c133,c211,c212,c213,c222,c223,c2&
     &33,c311,c312,c313,c322,c323,c333,k11,k12,k13,k21,k22,k23,k31,k32,k&
     &33,alpha,beta1,beta2,beta3,h,dt,nx,ny,nz)
      call evolve_k (gup11,gup12,gup13,gup22,gup23,gup33,k11,k12,k13,k21&
     &,k22,k23,k31,k32,k33,k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32&
     &m1,k33m1,c111,c112,c113,c122,c123,c133,c211,c212,c213,c222,c223,c2&
     &33,c311,c312,c313,c322,c323,c333,r11,r12,r13,r22,r23,r33,alpha,bet&
     &a1,beta2,beta3,nx,ny,nz,h,dt)
      dt = dtsave
      time = timesave
      return
      end subroutine

