diff --git a/NOBUILD/Dev/01.test_codes.R b/NOBUILD/Dev/01.test_codes.R
index 5288a26..0cd1a24 100644
--- a/NOBUILD/Dev/01.test_codes.R
+++ b/NOBUILD/Dev/01.test_codes.R
@@ -79,13 +79,16 @@ d. <- d.[d. < max(d.)]
 Du <- lapply(d., do_calc_by_depth)
 #str(Du,1)
 
-recomp <- FALSE
+recomp <- TRUE
 
+qt.nbgf <- resQt(m_dot=2e9/(5*365*86400), # 2014--2021
+		resRadius=hw, 
+		timespan=12*30*86400)
 if (!exists('cu.o') | recomp){
 cu.o <- pbapply::pbsapply(Du, cfs_isoporo, 
 	Beta.deg=opt_ang, fric=fric, Skemp=skemp, 
 	Shearmod = 10e9, nu_u = poiss, 
-	Qt = resQt(m_dot=1e6/(12*30*86400), resRadius=hw, timespan=12*30*86400),
+	Qt = qt.nbgf,
 	verbose=FALSE)
 }
 
@@ -94,7 +97,8 @@ r0 <- raster::raster(cu.o)
 rF <- raster::focal(r0, 
 	w=matrix(1, nrow=3, ncol=3), # equal-weighted average of all surrounding points
 	fun=mean, NAonly=TRUE, na.rm=TRUE) 
-cu <- as.matrix(rF)
+
+cu <- as.matrix(rF) / 1e6 # convert to ?Pa
 
 #cu <- sapply(Du, mean_stress)
 #cu <- sapply(Du, max_shear)
@@ -116,8 +120,9 @@ lCup <- c(xy, list(z = log10(cu)))
 lCun <- c(xy, list(z = log10(-cu)))
 
 fields::image.plot(Cu, asp=1, ylim=rev(range(Cu[['y']], na.rm=TRUE)),
-	col = c(NA,viridis::cividis(64, direction=-1)), 
-	zlim=max(abs(cu),na.rm=TRUE)*c(0,1))
+	col = c(NA,viridis::cividis(64, direction=-1)))
+	#, 
+	#zlim=max(abs(cu),na.rm=TRUE)*c(-1,1))
 #image(Cupos, add=TRUE, col=adjustcolor(c('white',NA), alpha=0.5))
 #contour(lCup, add=TRUE)
 #contour(lCun, lty=3, add=TRUE)
diff --git a/NOBUILD/Dev/LambertCodes.R b/NOBUILD/Dev/LambertCodes.R
index ea9719f..3db60ea 100644
--- a/NOBUILD/Dev/LambertCodes.R
+++ b/NOBUILD/Dev/LambertCodes.R
@@ -196,6 +196,7 @@ LambertTsai2020_Uniform <- function(x1, x2, resTopDepth, resThickness, resHalfWi
 	ret
 }
 
+
 #--------------------------------------------------------------
 # Plane-stress computations
 
@@ -249,6 +250,7 @@ planestress_principal_stresses.LT.diffusive <- planestress_principal_stresses.LT
 		vals
 	}
 	s1s2s3 <- t(apply(S, 1, .get_eig))
+	print(summary(s1s2s3))
 	colnames(s1s2s3) <- c('PS1','PS2','PS3')
 	cbind(s11=s1, s12=s12, s22=s2, s1s2s3)
 }
@@ -276,14 +278,20 @@ prefactor.LT.uniform <- function(x, Shearmod, nu_u, Skemp, rho_f, Qt=1, verbose=
 		rho_f <- 1000.0
 	}
 	if (!is.null(warns)) warning(warns, call.=FALSE)
-	# between eqs 8 / 9
+	#
+	# Scaling in Equation 14
+	#
+	resRadius <- x[['res']]$half.width
+	C <- Shearmod * (1 + nu_u) * Skemp / (2 * resRadius * 3 * pi * rho_f * (1 - nu_u))
+	
 	# result still needs to be multiplied by Q*t
 	# where Q is net mass flux from the producing region, t is time
 	# Q = dm/dt / Area = rho_f * volume
 	# dm/dt = rho_f * volume * area
 	# good simple description: https://web.mit.edu/16.unified/www/FALL/fluids/Lectures/f06.pdf
-	C <- Shearmod * (1 + nu_u) * Skemp / (3 * pi * rho_f * (1 - nu_u))
-	C * Qt / 2
+	
+	C * Qt
+	
 }
 prefactor.LT.diffusive <- function(x, ...) .NotYetImplemented()
 
@@ -311,11 +319,13 @@ fault_stresses.LT.diffusive <- fault_stresses.LT.uniform <- function(x, Beta=45,
 cfs_isoporo.LT.diffusive <- cfs_isoporo.LT.uniform <- function(x, Beta.deg, fric, Skemp, verbose=TRUE, ...){
 	if (verbose) message('CFS {S1-to-fault angle = ', Beta.deg, '°, friction = ', fric, "} for isotropic undrained poroelastic response {B = ", Skemp, "}")
 	nfs <- fault_stresses(x, Beta=Beta.deg, is.deg=TRUE)
+	#print(summary(nfs))
 	C <- prefactor(x, Skemp = Skemp, ...)
 	fs <- C * nfs
 	Tau <- fs[,'Tau']
 	Snorm <- fs[,'Snorm']
 	Smean <- fs[,'Smean']
-	cfsip <- cfs_isoporo.default(tau=Tau, mu=fric, Snormal = Snorm, Smean = Smean, B=Skemp)
+	cfsip <- cfs_isoporo.default(tau=Tau, mu=fric, Snormal = Snorm, Smean = Smean, B = Skemp)
+	#print(summary(cfsip))
 	zapsmall(cfsip)
 }
diff --git a/NOBUILD/Dev/beeler_cfs.R b/NOBUILD/Dev/beeler_cfs.R
index 79e738a..438d4e8 100644
--- a/NOBUILD/Dev/beeler_cfs.R
+++ b/NOBUILD/Dev/beeler_cfs.R
@@ -13,7 +13,10 @@ cfs_basic <- function(tau, mu, Snormal, pp){
 cfs_isoporo <- function(x, ...) UseMethod("cfs_isoporo")
 cfs_isoporo.default <- function(tau, mu, Snormal, Smean, B, ...){
 	pp.iso <- pp_isotropic(Smean, B)
-	cfs_basic(tau, mu, Snormal, pp.iso)
+	cf <- cfs_basic(tau, mu, Snormal, pp.iso)
+	print("")
+	print(summary(cbind(Smean, B, pp.iso, tau, cf)))	
+	cf
 }
 
 apparent_friction_basic <- function(mu, pp, Snormal){