      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


      subroutine constraint(nx, ny, nz, k11, k12, k13, k21, k22, k23, k3&
     &1, k32, k33, c111, c112, c113, c122, c123, c133, c211, c212, c213,&
     & c222, c223, c233, c311, c312, c313, c322, c323, c333, r11, r12, r&
     &13, r22, r23, r33, gup11, gup12, gup13, gup22, gup23, gup33, hamno&
     &rm, 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
      real (SP),  dimension(size (k11,1),size (k11,2),size (k11,3)):: t1&
     &,t2,t3,t4,t5,t6,t7,t8,t9,temp1 
!HPF$      align with k11 :: t1,t2,t3,t4,t5,t6,t7,t8,t9,temp1
      integer (I4B):: i,j
      t7 = k11+k22+k33
      temp1 = gup11*r11+gup12*r12*2.0d0+gup13*r13*2.0d0+gup22*r22+gup23*&
     &r23*2.0d0+gup33*r33+t7*t7+(k11*k11+k12*k21+k13*k31+k21*k12+k22*k22&
     &+k23*k32+k31*k13+k32*k23+k33*k33)
      call diff_x (nx,ny,nz,h,k11,t1)
      call diff_y (nx,ny,nz,h,k21,t2)
      call diff_z (nx,ny,nz,h,k31,t3)
      call diff_x (nx,ny,nz,h,k11,t4)
      call diff_x (nx,ny,nz,h,k22,t5)
      call diff_x (nx,ny,nz,h,k33,t6)
      t7 = ((t1+t2+t3+(c111*k11+c112*k21+c113*k31+c212*k11+c222*k21+c223&
     &*k31+c313*k11+c323*k21+c333*k31))-(c111*k11+c211*k12+c311*k13+c112&
     &*k21+c212*k22+c312*k23+c113*k31+c213*k32+c313*k33))-(t4+t5+t6)
      call diff_x (nx,ny,nz,h,k12,t1)
      call diff_y (nx,ny,nz,h,k22,t2)
      call diff_z (nx,ny,nz,h,k32,t3)
      call diff_y (nx,ny,nz,h,k11,t4)
      call diff_y (nx,ny,nz,h,k22,t5)
      call diff_y (nx,ny,nz,h,k33,t6)
      t8 = ((t1+t2+t3+(c111*k12+c112*k22+c113*k32+c212*k12+c222*k22+c223&
     &*k32+c313*k12+c323*k22+c333*k32))-(c112*k11+c212*k12+c312*k13+c122&
     &*k21+c222*k22+c322*k23+c123*k31+c223*k32+c323*k33))-(t4+t5+t6)
      call diff_x (nx,ny,nz,h,k13,t1)
      call diff_y (nx,ny,nz,h,k23,t2)
      call diff_z (nx,ny,nz,h,k33,t3)
      call diff_z (nx,ny,nz,h,k11,t4)
      call diff_z (nx,ny,nz,h,k22,t5)
      call diff_z (nx,ny,nz,h,k33,t6)
      t9 = ((t1+t2+t3+(c111*k13+c112*k23+c113*k33+c212*k13+c222*k23+c223&
     &*k33+c313*k13+c323*k23+c333*k33))-(c113*k11+c213*k12+c313*k13+c123&
     &*k21+c223*k22+c323*k23+c133*k31+c233*k32+c333*k33))-(t4+t5+t6)
      hamnorm = maxval (temp1)
      mom1norm = maxval (t7)
      mom2norm = maxval (t8)
      mom3norm = maxval (t9)
      return
      end subroutine

