Viscosity

Viscosity is included mainly for numerical analysis so I chose $\sigma$ to be a constant. To create viscosity, the heat flow equation can be applied to the pressure. Unfortunately we cannot use the heat flow subroutine itself because the viscosity subroutine also needs to pass through the velocity. Starting from the heatflow subroutine, including the pass through, and increasing from one to two dimensions gives subroutine viscosity()

# Viscosity operator = heatflow on pressure and pass through velocity
#
subroutine viscosity( adj, add, sigma, qq,nx,ny,  rr   )
integer               adj, add,           nx,ny,        x, y,     p,   u,   v
real		 	        sigma, qq(nx,ny,3), rr(nx,ny,3)
							       p=1; u=2; v=3
call        adjnull(  adj, add,        qq,nx*ny*3,  rr,nx*ny*3)
do x= 2, nx-1 {
do y= 2, ny-1 {
	if( adj == 0) {
		rr(x,y,p) = rr(x,y,p) + qq(x,y,p) +
			( qq(x-1,y  ,p) - 2*qq(x,y,p) + qq(x+1,y  ,p)) * sigma +
			( qq(x  ,y-1,p) - 2*qq(x,y,p) + qq(x  ,y+1,p)) * sigma
		rr(x,y,u) = rr(x,y,u) + qq(x,y,u)
		rr(x,y,v) = rr(x,y,v) + qq(x,y,v)
	} else {
		qq(x,y,u) = qq(x,y,u) + rr(x,y,u)
		qq(x,y,v) = qq(x,y,v) + rr(x,y,v)
		qq(x,y,p) = qq(x,y,p) + rr(x,y,p) * ( 1 - sigma * 4)
		qq(x-1,y,p) = qq(x-1,y,p) + rr(x,y,p) * sigma
		qq(x+1,y,p) = qq(x+1,y,p) + rr(x,y,p) * sigma
		qq(x,y-1,p) = qq(x,y-1,p) + rr(x,y,p) * sigma
		qq(x,y+1,p) = qq(x,y+1,p) + rr(x,y,p) * sigma
		}
	}}
return; end