
subroutine nusselt(th,  n, m, condnf, condf, avnb, fii , dp, x1, x2, x3, y1, y2)
real th(0:n,0:m)
! Nusselt number Calculation
snub=0.0
do i=x1,x2
curth=th(i,0)
call condcalc(curth,condrate, fii, dp)
rnub=-(th(i,1)-th(i,0))*float(n)*condrate
snub=snub+rnub
enddo

do i=x3,n
curth=th(i,0)
call condcalc(curth,condrate, fii, dp)
rnub=-(th(i,1)-th(i,0))*float(n)*condrate
snub=snub+rnub
enddo

do j=y1,y2
curth=th(n,j)
call condcalc(curth,condrate, fii, dp)
rnub=-(th(n-1,j)-th(n,j))*float(n)*condrate
snub=snub+rnub
enddo

avnb=snub/(x2-x1+n-x3+y2-y1)


return
end