      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 riccidd(r11, r12, r13, r22, r23, r33, c111, c112, c113,&
     & c122, c123, c133, c211, c212, c213, c222, c223, c233, c311, c312,&
     & c313, c322, c323, c333, nx, ny, nz, h)
      use types
      use hpfprocessors
      use diff_interface
      real (SP),  dimension(:,:,:),intent(out):: r11,r12,r13,r22,r23,r33&
     & 
      real (SP),  dimension(:,:,:),intent(in):: c111,c112,c113,c122,c123&
     &,c133 
      real (SP),  dimension(:,:,:),intent(in):: c211,c212,c213,c222,c223&
     &,c233 
      real (SP),  dimension(:,:,:),intent(in):: c311,c312,c313,c322,c323&
     &,c333 
!HPF$      distribute r11 * (*,*,block) onto * PROC
!HPF$      align with * r11 :: r12,r13,r22,r23,r33
!HPF$      align with * r11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with * r11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with * r11 :: c311,c312,c313,c322,c323,c333
      integer (I4B):: nx,ny,nz
      real (SP):: h
      real (SP),  dimension(size (r11,1),size (r11,2),size (r11,3)):: t1&
     &1,t12,t13,t22,t23,t33 
!HPF$      align with r11 :: t11,t12,t13,t22,t23,t33
      integer (I4B):: i
      call diff_x (nx,ny,nz,h,c111,r11)
      call diff_x (nx,ny,nz,h,c112,r12)
      call diff_x (nx,ny,nz,h,c113,r13)
      call diff_x (nx,ny,nz,h,c122,r22)
      call diff_x (nx,ny,nz,h,c123,r23)
      call diff_x (nx,ny,nz,h,c133,r33)
      call diff_y (nx,ny,nz,h,c211,t11)
      call diff_y (nx,ny,nz,h,c212,t12)
      call diff_y (nx,ny,nz,h,c213,t13)
      call diff_y (nx,ny,nz,h,c222,t22)
      call diff_y (nx,ny,nz,h,c223,t23)
      call diff_y (nx,ny,nz,h,c233,t33)
      r11 = r11+t11
      r12 = r12+t12
      r13 = r13+t13
      r22 = r22+t22
      r23 = r23+t23
      r33 = r33+t33
      call diff_z (nx,ny,nz,h,c311,t11)
      call diff_z (nx,ny,nz,h,c312,t12)
      call diff_z (nx,ny,nz,h,c313,t13)
      call diff_z (nx,ny,nz,h,c322,t22)
      call diff_z (nx,ny,nz,h,c323,t23)
      call diff_z (nx,ny,nz,h,c333,t33)
      r11 = r11+t11
      r12 = r12+t12
      r13 = r13+t13
      r22 = r22+t22
      r23 = r23+t23
      r33 = r33+t33
      call diff_x (nx,ny,nz,h,c111,t11)
      call diff_y (nx,ny,nz,h,c111,t12)
      call diff_z (nx,ny,nz,h,c111,t13)
      call diff_y (nx,ny,nz,h,c112,t22)
      call diff_z (nx,ny,nz,h,c112,t23)
      call diff_z (nx,ny,nz,h,c113,t33)
      r11 = r11-t11
      r12 = r12-t12
      r13 = r13-t13
      r22 = r22-t22
      r23 = r23-t23
      r33 = r33-t33
      call diff_x (nx,ny,nz,h,c212,t11)
      call diff_y (nx,ny,nz,h,c212,t12)
      call diff_z (nx,ny,nz,h,c212,t13)
      call diff_y (nx,ny,nz,h,c222,t22)
      call diff_z (nx,ny,nz,h,c222,t23)
      call diff_z (nx,ny,nz,h,c223,t33)
      r11 = r11-t11
      r12 = r12-t12
      r13 = r13-t13
      r22 = r22-t22
      r23 = r23-t23
      r33 = r33-t33
      call diff_x (nx,ny,nz,h,c313,t11)
      call diff_y (nx,ny,nz,h,c313,t12)
      call diff_z (nx,ny,nz,h,c313,t13)
      call diff_y (nx,ny,nz,h,c323,t22)
      call diff_z (nx,ny,nz,h,c323,t23)
      call diff_z (nx,ny,nz,h,c333,t33)
      r11 = r11-t11
      r12 = r12-t12
      r13 = r13-t13
      r22 = r22-t22
      r23 = r23-t23
      r33 = r33-t33
      r11 = r11+(c111+c212+c313)*c111+(c112+c222+c323)*c211+(c113+c223+c&
     &333)*c311
      r12 = r12+(c111+c212+c313)*c112+(c112+c222+c323)*c212+(c113+c223+c&
     &333)*c312
      r13 = r13+(c111+c212+c313)*c113+(c112+c222+c323)*c213+(c113+c223+c&
     &333)*c313
      r22 = r22+(c111+c212+c313)*c122+(c112+c222+c323)*c222+(c113+c223+c&
     &333)*c322
      r23 = r23+(c111+c212+c313)*c123+(c112+c222+c323)*c223+(c113+c223+c&
     &333)*c323
      r33 = r33+(c111+c212+c313)*c133+(c112+c222+c323)*c233+(c113+c223+c&
     &333)*c333
      r11 = ((r11-(c111*c111+c211*c112+c311*c113))-(c112*c211+c212*c212+&
     &c312*c213))-(c113*c311+c213*c312+c313*c313)
      r12 = ((r12-(c111*c112+c211*c122+c311*c123))-(c112*c212+c212*c222+&
     &c312*c223))-(c113*c312+c213*c322+c313*c323)
      r13 = ((r13-(c111*c113+c211*c123+c311*c133))-(c112*c213+c212*c223+&
     &c312*c233))-(c113*c313+c213*c323+c313*c333)
      r22 = ((r22-(c112*c112+c212*c122+c312*c123))-(c122*c212+c222*c222+&
     &c322*c223))-(c123*c312+c223*c322+c323*c323)
      r23 = ((r23-(c112*c113+c212*c123+c312*c133))-(c122*c213+c222*c223+&
     &c322*c233))-(c123*c313+c223*c323+c323*c333)
      r33 = ((r33-(c113*c113+c213*c123+c313*c133))-(c123*c213+c223*c223+&
     &c323*c233))-(c133*c313+c233*c323+c333*c333)
      end subroutine

