diff --git a/src/bk.f b/src/bk.f
index e2f0faf0595e31f08b4417036bc0c72e56bf5eb3..25f7eeb7246b3d8d3d00bff234574282f10dc2f9 100644
--- a/src/bk.f
+++ b/src/bk.f
@@ -252,7 +252,7 @@ c     determine search neighbourhood
       CALL SRCHNB(LON0,LAT0,DO0,INDDO,N0,NDO,LON,LAT,N,NS,
      .            RSEARCH,NSEARCH,NSMIN,NSMAX,
      .            INDSRT,INDSNB,INDSNW,DIST)
-      write(*,*)"SEARCHNB is:",NS
+c      write(*,*)"SEARCHNB is:",NS
       if (usesbbt.ne.0) then
          do 1212 i=1,n
             if (indsnb(i).ne.0) then
@@ -416,9 +416,9 @@ c     prepare the right hand side(s)
       
 c     solve the system
 c     LAMBDA = (C0VEC+FWORK*PHIPR*F0WORK)*(COVMAT+FWORK*PHIPR*FWORK')^-1
-      name="kwork\0"
-      call matpr(name,kwork,ns,ns,LDkwrk,dbglvl)
-      name="rhswork\0"
+c      name="kwork\0"
+c      call matpr(name,kwork,ns,ns,LDkwrk,dbglvl)
+c      name="rhswork\0"
       call matpr(name,rhswork,ns,NDO,LDkWrk,dbglvl)
 C     BETTER USE DSYSVX ?
 
@@ -489,15 +489,17 @@ c                                                                    C0VEC)
       IF (MODE.EQ.2 .OR. MODE.EQ.3) THEN
 c            work array: RHSWORK = C0VEC + FWORK PHIPR F0WORK
 c     the same procedure as to prepare the right hand side(s) above:
+c     (use the same index sorting as in lambda!!)
       DO 1300 J=1,NDO
          DO 140 I=1,NS
-            RHSWORK(I,J)=C0VEC(INDSNB(I),inddo(J))+FPF0WRK(I,j)
+            RHSWORK(indsnb(I),J)=C0VEC(INDSNB(I),inddo(J))+FPF0WRK(I,J)
  140     CONTINUE
  1300 CONTINUE
          DO 500  I=1,NDO
 c           work array: MU = PHIPR * F0WORK
             CALL DGEMV('N',NTREND,NTREND,ONE,PHIWRK,LDPHWK,
-     .                 F0WORK(1,I),1,ZERO,MU,1)         
+     .                 F0WORK(1,I),1,ZERO,MU,1)        
+c         write(*,*)"cov0:",cov0
 c         name="mu\0"
 c         call matpr(name,mu,ntrend,1,ldmu,dbglvl)
 c         name="f0work\0"
@@ -508,6 +510,8 @@ c         name="lambda\0"
 c         call matpr(name,lambda(1,inddo(i)),n,1,ldlmbd,dbglvl)
 c         name="c0vec\0"
 c         call matpr(name,c0vec(1,inddo(i)),n,1,ldc0,dbglvl)
+c         name="rhswork\0"
+c         call matpr(name,rhswork(1,inddo(i)),n,1,ldkwrk,dbglvl)
               
             VAR0(INDDO(I)) = COV0 - DDOT(N,LAMBDA(1,INDDO(I)),
      .                                   1,RHSWORK(1,INDDO(I)),1)
diff --git a/src/bk_point.c b/src/bk_point.c
index 3f65a1dbf167fa4749bac4f20d5b37830c4e4763..38c7409adfd4a5a1e397c452211170f451857332 100644
--- a/src/bk_point.c
+++ b/src/bk_point.c
@@ -137,9 +137,9 @@ void bk_point(double *xp,
 		 dop,
 		 inddop,
 		 np,
-		     lon,
-		     lat,
-		     z,
+		 lon,
+		 lat,
+		 z,
 		     n,
 		     covtype,
 		     covpar,